Method and device for simulating freezing process of saline water film on surface of typical component of polar ship

Abstract

The invention discloses a method and device for simulating the freezing process of a salt water film on the surface of a typical component of a polar ship, belongs to the technical field of salt water freezing process simulation, and particularly relates to simulation of the freezing process of the salt water film on the surface of the typical component of the polar ship. The problems that research on the existing brine freezing process is too simplified, and factors such as dynamic influence of salt migration (precipitation) on the freezing process and dynamic change of a contact angle caused by a hysteresis effect are not considered, so that simulation of the seawater droplet freezing process is not accurate enough are solved. The method comprises a solving step used for solving a momentum equation, an energy equation and a concentration equation according to input parameters of a saline film freezing numerical model to obtain output parameters of the saline film freezing numerical model. The method and the device for simulating the freezing process of the saline water film on the surface of the typical component of the polar ship are suitable for simulating the freezing process of the saline water film on the surface of the typical component of the polar ship.

Description

Technical Field

[0001] The invention relates to the technical field of simulation of salt water freezing process, and in particular to simulation of the freezing process of salt water film on the surface of typical components of polar ships. Background Art

[0002] Ships sailing in high latitudes such as the Arctic often face icing problems on their superstructures. Large amounts of ice accumulation can affect ship stability, hinder the normal operation of equipment, and even endanger personnel safety. Taking the icing of polar ships as an example, the process is as follows: Figure 4 As shown, a large number of droplets will be formed after the bow collides with the waves. Some of the droplets move under their own gravity and external forces and cover the superstructure. Due to the presence of salt in seawater, the supercooled salt water droplets in the droplets will only freeze into solid particles in the air when the temperature is lower than -23°C, but the air temperature is usually higher than -23°C, so the water droplets in the droplets are mostly liquid in the air. When the air temperature is lower than the freezing temperature of seawater (about -2°C) and the supercooled salt water droplets hit the surface of the structure, the droplets will begin to freeze. The cold air carries a large number of liquid water droplets (droplets) and hits the surface of the structure to form a liquid water film. Most of the water in the water film will leave the surface in the form of "run-off" under the action of gravity, and a small amount of water will adhere to the surface of the structure to form ice accumulation. This ice growth process occurring on the surface of the structure is usually called "wet growth".

[0003] The freezing of seawater is accompanied by the growth of pure water ice crystals, the migration and precipitation of salt. The salt in the gaps between ice crystals gathers to form brine cells, and multiple brine cells are connected to form brine channels (high-salinity solution channels). The brine channels are irregularly distributed along the gaps between ice crystals, which will affect the strength of sea ice. When the volume fraction of brine increases from 10% to 30%, the corresponding shear strength and bending strength of sea ice decrease by 70% and 64% respectively.

[0004] The freezing of seawater droplets is a complex thermophysical process that includes macroscopic phenomena (scale greater than 1 mm) such as brine film runoff, ice layer growth, and flow of high-concentration salt solutions, as well as microscopic phenomena (scale between 1 nm and 1 mm) such as supercooled water droplet impact freezing, ice crystal nucleation growth, and salt precipitation. It is also affected by factors such as the inclination angle and wettability of the surface of the ship's superstructure structure. Its accurate simulation requires overcoming technical difficulties such as salt removal modeling, gas-liquid-solid coupling, and analysis of supercooled water droplet dynamics and thermodynamic behavior.

[0005] Existing studies on salt water freezing (or seawater droplet freezing models) have the following problems: (1) Most existing seawater droplet freezing models are based on empirical formulas, theories or commercial software, focusing on the prediction of macroscopic large-scale ice thickness, but not on important phenomena that affect the final ice morphology, such as the evolution of brine channels and runoff inside the brine film during the freezing process. This results in rough and qualitative prediction results. (2) Existing studies on water films at the mesoscopic scale are mostly focused on horizontal surface scenarios. However, the surface of ship superstructure structures also contains a large number of inclined and vertical structural surfaces. The study of brine films on horizontal surfaces alone cannot reflect the factors of brine channel orientation under the action of runoff and gravity. Therefore, it is necessary to develop a method that can simulate the freezing of brine films on isolated cold plates with angles, so as to improve the research system of ship superstructure icing prediction. (3) Most existing studies on salt water freezing only consider the static effect of salinity on the freezing point. That is, the freezing temperature of pure water is 0°C. After considering salinity, the corresponding freezing temperature is lower than 0°C, which in turn affects the freezing process. This setting can only reflect the static effect of salinity, while ignoring or simplifying the dynamic effect of salt migration on the freezing process, which limits the refined prediction and evaluation of seawater droplet freezing. (4) The water droplets in the polar ship icing environment are supercooled salt water droplets (i.e., they can remain in the liquid phase when the temperature is lower than the freezing temperature). Supercooled salt water droplets in metastable equilibrium only need a small disturbance to trigger solidification and become steady-state solid ice. The difference between the temperature of the supercooled salt water droplets and the freezing point is called supercooling. Supercooling affects processes such as ice crystal nucleation when the water droplets freeze. Moreover, when a single supercooled water droplet moves to the surface of the structure, it may exhibit a variety of dynamic and thermodynamic behaviors such as no collision, collision and spreading, retraction and bounce, bottom sliding, overflow, adhesion and freezing under different surface wettability, air temperature and movement speed. Therefore, factors such as water droplet supercooling, surface tension, heat and initial kinetic energy need to be considered in numerical simulation. However, in the traditional model, the water droplet is simplified into a spherical unit with constant shape and mass, and these factors are not considered. (5) The different surface roughness of different materials will cause surface hysteresis effect. In the simulation of water droplet impact on the surface of the structure, it is also necessary to consider the dynamic change of contact angle caused by the hysteresis effect. These problems are not considered in the existing model.

[0006] In summary, the existing research on the brine freezing process is overly simplified. It lacks consideration of factors such as the evolution and runoff of brine channels inside the brine film during the freezing process, inclined or vertical structural surfaces, water droplet supercooling, surface tension, heat and initial kinetic energy, the dynamic impact of salt migration on the freezing process, and the change in dynamic contact angle caused by the hysteresis effect. This leads to the problem of inaccurate simulation of the seawater droplet freezing process. Summary of the invention

[0007] The present invention proposes a method and device for simulating the freezing process of brine film on the surface of typical components of polar ships, which solves the problem that the existing research on the freezing process of brine is overly simplified, and lacks consideration of factors such as the evolution and runoff of brine channels inside the brine film during the freezing process, inclined or vertical structural surfaces, water droplet supercooling, surface tension, heat and initial kinetic energy, the dynamic influence of salt migration on the freezing process, and the change of dynamic contact angle caused by the hysteresis effect, resulting in insufficient accuracy in the simulation of the freezing process of seawater droplets.

[0008] The method for simulating the freezing process of salt water film on the surface of typical components of polar ships described in the present invention has the following technical solutions: The method comprises the following steps: Input step: used to obtain input parameters of the salt water film icing numerical model according to the characteristics of typical components of polar navigation ships and navigation environment information; The input parameters of the numerical model of brine film icing include physical parameters and calculation parameters; The physical property parameters include seawater density and salinity, liquid water content, dynamic viscosity, initial temperature, eutectic concentration, eutectic temperature, partition coefficient, liquid phase latent heat, specific heat capacity, thermal diffusion coefficient, component diffusion coefficient, expansion coefficient, back diffusion coefficient, pastiness coefficient, latent heat renewal relaxation factor and liquid phase fraction; The calculation parameters include calculation area, number of grids, time step, initial conditions, boundary conditions, water drop release interval and number of water drops; Solving step: solving the momentum equation, energy equation and concentration equation according to the input parameters of the brine film icing numerical model, and obtaining the output parameters of the brine film icing numerical model; The momentum equation is a fluid equation established in accordance with Newton's second law; The energy equation is an energy conservation equation, which complies with the first law of thermodynamics and takes temperature as the variable to be determined; The concentration equation is a component conservation equation that complies with the component mass conservation law and is used to obtain the concentration parameter value in the calculation domain; The output parameters are temperature, salinity, velocity, pressure, liquid fraction, volume fraction, signed distance function and mixing fraction; The solving steps include a dynamic contact angle application step, an isolated cold plate step, a solid-liquid interface tracking step, a gas-liquid interface tracking step, a microscopic anti-diffusion model selection step, and a general component step; The dynamic contact angle application step is used to set the contact angle boundary conditions on the cold plate surface to reflect the different wettabilities on the cold plate surface and to consider the surface hysteresis effect caused by the different wettabilities; The isolated cold plate step is used to set an isolated cold plate area that is not connected to the outer boundary in the air calculation domain to simulate the runoff phenomenon of the salt water film; The solid-liquid interface tracking step is used to solve the energy equation and obtain the temperature and liquid fraction of each grid node in the calculation area; The gas-liquid interface tracking step is used to solve the volume fraction and signed distance function that characterize the gas-liquid interface, and select the VOSET gas-liquid interface tracking method to process the gas-liquid interface changes during the motion of the water droplet; it is also used to add a source term about the liquid phase fraction to the momentum equation, and by setting the coefficient of the source term, the fluidity is maintained when the liquid phase ratio is high, and the flow is suppressed when the liquid phase ratio is low; it is also used to solve the momentum equation using the IDEAL algorithm based on the velocity pressure correction theory; The microscopic anti-diffusion model selection step is used to solve the concentration equation using a fully implicit time integration scheme; and is also used to select an anti-diffusion model to characterize the concentration distribution; The general component step includes an algebraic equation group step, a control equation discrete coefficient step, a grid coordinate step, and a coefficient zeroing step, wherein: The algebraic equation group step is used to solve the discrete momentum, temperature and concentration control equations, and is used to be called by the solid-liquid interface tracking step and the microscopic anti-diffusion model selection step; The control equation discrete coefficient step is used to solve the coefficients of the velocity and pressure equations in the momentum equation, thereby providing input for the algebraic equation system step; The grid coordinate step is used to divide the calculation area and define the grid coordinate values ​​in the calculation area; The coefficient zeroing step is used to zero the coefficients of the control equation before the control equation discrete coefficient step, so that the control equation discrete coefficient step can be repeatedly executed; Output step: used to obtain output results according to the output parameters of the brine film icing numerical model; the output results include the growth rate of the surface ice layer, ice accumulation characteristics and the influence of different factors on the ice accumulation characteristics.

[0009] Furthermore, a preferred embodiment is provided, in which the value of the mixture fraction in the output parameter is half of the sum of the liquid phase fraction and the volume fraction, which is used to solve the problem that the liquid phase fraction or the volume fraction alone cannot simultaneously reflect and distinguish the three phases of solid, liquid and gas.

[0010] Further, a preferred embodiment is provided, wherein the setting of the cold plate surface contact angle boundary condition to reflect the different wettabilities of the cold plate surface and the surface hysteresis effect caused by the different wettabilities includes: The dynamic contact angle of the water droplet during its motion is solved in real time, and the dynamic contact angle input is provided for the gas-liquid interface tracking step, thereby correcting the signed distance function value at the solid boundary, thereby reflecting the different wettability of the cold plate surface in the numerical model of brine film icing; The different wettabilities include hydrophilicity, hydrophobicity and superhydrophobicity; The hydrophilicity is a contact angle of less than 90°; The hydrophobicity is a contact angle between 90° and 150°; The super hydrophobicity is a contact angle greater than 150°; The surface hysteresis effect caused by different wettabilities is the dynamic change of the contact angle between the water droplet and the solid surface during the spreading and retracting process.

[0011] Further, a preferred embodiment is provided, wherein the real-time solution of the dynamic contact angle during the motion of the water droplet comprises: The Kistler dynamic contact angle model is used to solve the dynamic contact angle. :

[0012] In the formula, is the Hoffman smoothing function, which is calculated as follows:

[0013] In the formula, the independent variable Defined as:

[0014] In the formula, is the inverse function of the Hoffman smoothing function, which is solved by the following iteration:

[0015] in, is the capillary number, which is used to characterize the ratio of viscosity to surface tension; To balance the contact angle, The inverse function of , used to reflect the contact angle hysteresis effect; According to the three-dimensional three-phase contact point velocity of water droplet-air-fixed surface direction to determine :

[0016] In the actual scene of water drop movement, the following formula is used to solve :

[0017] in,

[0018] in, is the distance between the three-phase contact point and the initial impact point of the water droplet.

[0019] Further, a preferred embodiment is provided, in the step of isolating the cold plate: The isolated cold plate is included in the fluid calculation domain and solved as a whole, so that the velocity in the isolated cold plate area is zero and the temperature is a given sub-zero value. : The large coefficient method is used to make the temperature in the isolated cold plate area take a given value : Before solving the temperature algebraic equations, let the nodes in the isolated cold plate region The corresponding discrete equation is: ; In the formula, , ,in, is not less than a larger value of ; To isolate the nodes in the cold plate area Corresponding temperature Always take the given value ; Will The corresponding position is set to 0, so that the velocity field in the isolated cold plate area is always zero; To avoid pressure correction to velocity in isolated cold plate regions, which would destroy the zero velocity condition, the following settings are repeated in each time step: The first setting is: before each iteration, the velocity in the isolated cold plate area is set to 0; because of the staggered grid, when When is the boundary node in the isolated cold plate area, let , similarly When is the right boundary, ; The second setting: Before solving the velocity algebraic equations, let the main diagonal coefficients of each velocity discrete matrix in the isolated cold plate region be is not less than A larger value of , thereby ensuring the calculated intermediate speed ; The third setting: Use the IDEAL algorithm to solve the pressure-velocity coupling. Calculate the pressure correction value in two inner iterations. To avoid the pressure correction value in the isolated cold plate area to correct the velocity, set the velocity correction equation coefficient to be not greater than The minimum value of .

[0020] Further, a preferred embodiment is provided, in the solid-liquid interface tracking step, the temperature and liquid fraction of each grid node in the calculation area are obtained by the following method: The source term of the energy equation contains a transient term of liquid fraction; Iteratively solve the liquid fraction and obtain the temperature after the convergence result is obtained; When solving iteratively, only one iteration is required at a time: After the solution of the layer velocity and pressure at the current time is completed, the liquid phase fraction is updated and its upper and lower limits are specified, thereby obtaining the source term of the energy equation and the coefficient of the discrete equation; Using the steps of algebraic equations, solve the energy equation to obtain the temperature value, determine the source term and the discrete equation coefficient, and solve the algebraic equations to obtain the current temperature; The above steps are iterated several times until the temperature and liquid fraction of the current time layer reach consistent convergence conditions.

[0021] The present invention also proposes a device for simulating the freezing process of salt water film on the surface of typical components of polar ships, and its technical scheme is as follows: The device comprises the following modules: Input module: used to obtain the input parameters of the salt water film icing numerical model according to the characteristics of typical components of polar ships and navigation environment information; The input parameters of the numerical model of brine film icing include physical parameters and calculation parameters; The physical property parameters include seawater density and salinity, liquid water content, dynamic viscosity, initial temperature, eutectic concentration, eutectic temperature, partition coefficient, liquid phase latent heat, specific heat capacity, thermal diffusion coefficient, component diffusion coefficient, expansion coefficient, back diffusion coefficient, pastiness coefficient, latent heat renewal relaxation factor and liquid phase fraction; The calculation parameters include calculation area, number of grids, time step, initial conditions, boundary conditions, water drop release interval and number of water drops; Solver: used to solve the momentum equation, energy equation and concentration equation according to the input parameters of the brine film icing numerical model, and obtain the output parameters of the brine film icing numerical model; The momentum equation is a fluid equation established in accordance with Newton's second law; The energy equation is an energy conservation equation, which complies with the first law of thermodynamics and takes temperature as the variable to be determined; The concentration equation is a component conservation equation that complies with the component mass conservation law and is used to obtain the concentration parameter value in the calculation domain; The output parameters are temperature, salinity, velocity, pressure, liquid fraction, volume fraction, signed distance function and mixing fraction; The solution module includes a dynamic contact angle application component, an isolated cold plate component, a solid-liquid interface tracking component, a gas-liquid interface tracking component, a microscopic anti-diffusion model selection component and a general component; The dynamic contact angle applying component is used to set the contact angle boundary condition of the cold plate surface to reflect the different wettability of the cold plate surface and consider the surface hysteresis effect caused by the different wettability; The isolated cold plate component is used to set an isolated cold plate area that is not connected to the outer boundary in the air calculation domain to simulate the runoff phenomenon of the salt water film; The solid-liquid interface tracking component is used to solve the energy equation and obtain the temperature and liquid fraction of each grid node in the calculation area; The gas-liquid interface tracking component is used to solve the volume fraction and signed distance function that characterize the gas-liquid interface, and selects the VOSET gas-liquid interface tracking method to process the gas-liquid interface changes during the motion of the water droplet; it is also used to add a source term about the liquid phase fraction to the momentum equation, and by setting the coefficient of the source term, the fluidity is maintained when the liquid phase ratio is high, and the flow is suppressed when the liquid phase ratio is low; it is also used to solve the momentum equation using the IDEAL algorithm based on the velocity pressure correction theory; The microscopic anti-diffusion model selection component is used to solve the concentration equation using a fully implicit time integration scheme; it is also used to select an anti-diffusion model to characterize the concentration distribution; The general component includes an algebraic equation group module, a control equation discrete coefficient module, a grid coordinate module and a coefficient zeroing module, wherein: The algebraic equation group module is used to solve the discrete momentum, temperature and concentration control equations, and is used to be called by the solid-liquid interface tracking step and the microscopic anti-diffusion model selection step; The control equation discrete coefficient module is used to solve the coefficients of the velocity and pressure equations in the momentum equation, thereby providing input for the algebraic equation system step; The grid coordinate module is used to divide the calculation area and define the grid coordinate values ​​in the calculation area; The coefficient zeroing module is used to zero the coefficients of the control equation before the discrete coefficient step of the control equation, so that the discrete coefficient step of the control equation can be repeatedly executed; Output module: used to obtain output results according to the output parameters of the brine film icing numerical model; the output results include the growth rate of the surface ice layer, ice accumulation characteristics and the influence of different factors on the ice accumulation characteristics.

[0022] The present invention also proposes a computer device, and its technical solution is as follows: A computer device comprises: a processor and a memory, wherein the memory is used to store executable instructions of the processor, and the processor is configured to execute the above-mentioned method for simulating the icing process of salt water film on the surface of typical components of polar ships by executing the executable instructions.

[0023] The present invention also proposes a computer storage medium, and its technical solution is as follows: A computer storage medium stores a computer program, and when the computer program is run, the method for simulating the freezing process of salt water film on the surface of a typical component of a polar ship is executed.

[0024] The present invention also proposes a computer program product, and its technical solution is as follows: A computer program product includes a computer program / instruction, which, when executed by a processor, implements the steps of the above-mentioned method for simulating the freezing process of salt water film on the surface of typical components of polar ships.

[0025] The present invention has the following beneficial effects: The method for simulating the freezing process of brine film on the surface of typical components of polar ships described in the present invention starts from the movement-freezing scenario of supercooled water droplets and, based on the "island" theory, establishes a numerical method for the freezing of brine film taking into account the supercooling effect and the hysteresis effect. The method can simulate the process of forming a brine film after multiple water droplets hit typical components, the freezing of pure water in the film, the discharge of salt, etc., revealing the formation and evolution mechanism of brine channels, providing a reference for improving the traditional numerical model of seawater droplet freezing, and improving the accuracy of the seawater droplet freezing forecasting method.

[0026] The method and device for simulating the freezing process of salt water film on the surface of typical components of polar ships described in the present invention are suitable for simulating the freezing process of salt water film on the surface of typical components of polar ships. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0028] Figure 1 It is a schematic structural diagram of a device for simulating the freezing process of salt water film on the surface of a typical component of a polar ship in one embodiment of the present invention; Figure 2 A schematic diagram of the iterative solution process of the dynamic contact angle in one embodiment of the present invention; Figure 3 In one embodiment of the present invention, the signed distance function within the grid at the boundary is Solution diagram; Figure 4 The diagram is a schematic diagram of the movement of droplet clusters and the freezing phenomenon of droplets in one embodiment of the present invention. DETAILED DESCRIPTION

[0029] In order to make the technical solutions and advantages of the present invention more clearly described, the specific implementation methods of the present invention will be further described in detail and completely in conjunction with the accompanying drawings. The various implementation methods described below are only part of the preferred solutions of the present invention, rather than all implementation methods; the various implementation methods described below are intended to explain the present invention and cannot be understood as limitations on the present invention; the reasonable combination of technical features defined in the various implementation methods of the present invention, and all other implementation methods obtained by ordinary technicians in the field without creative work based on the implementation methods of the present invention, are within the scope of protection of the present invention.

[0030] One implementation method provides a method for simulating the freezing process of salt water film on the surface of typical components of polar ships: The method comprises the following steps: Input step: used to obtain input parameters of the salt water film icing numerical model according to the characteristics of typical components of polar navigation ships and navigation environment information; The input parameters of the numerical model of brine film icing include physical parameters and calculation parameters; The physical property parameters include seawater density and salinity, liquid water content, dynamic viscosity, initial temperature, eutectic concentration, eutectic temperature, partition coefficient, liquid phase latent heat, specific heat capacity, thermal diffusion coefficient, component diffusion coefficient, expansion coefficient, back diffusion coefficient, pastiness coefficient, latent heat renewal relaxation factor and liquid phase fraction; The calculation parameters include calculation area, number of grids, time step, initial conditions, boundary conditions, water drop release interval and number of water drops; Solving step: solving the momentum equation, energy equation and concentration equation according to the input parameters of the brine film icing numerical model, and obtaining the output parameters of the brine film icing numerical model; The momentum equation is a fluid equation established in accordance with Newton's second law; The energy equation is an energy conservation equation, which complies with the first law of thermodynamics and takes temperature as the variable to be determined; The concentration equation is a component conservation equation that complies with the component mass conservation law and is used to obtain the concentration parameter value in the calculation domain; The output parameters are temperature, salinity, velocity, pressure, liquid fraction, volume fraction, signed distance function and mixing fraction; The solving steps include a dynamic contact angle application step, an isolated cold plate step, a solid-liquid interface tracking step, a gas-liquid interface tracking step, a microscopic anti-diffusion model selection step, and a general component step; The dynamic contact angle application step is used to set the contact angle boundary conditions on the cold plate surface to reflect the different wettabilities on the cold plate surface and to consider the surface hysteresis effect caused by the different wettabilities; The isolated cold plate step is used to set an isolated cold plate area that is not connected to the outer boundary in the air calculation domain to simulate the runoff phenomenon of the salt water film; The solid-liquid interface tracking step is used to solve the energy equation and obtain the temperature and liquid fraction of each grid node in the calculation area; The gas-liquid interface tracking step is used to solve the volume fraction and signed distance function that characterize the gas-liquid interface, and select the VOSET gas-liquid interface tracking method to process the gas-liquid interface changes during the motion of the water droplet; it is also used to add a source term about the liquid phase fraction to the momentum equation, and by setting the coefficient of the source term, the fluidity is maintained when the liquid phase ratio is high, and the flow is suppressed when the liquid phase ratio is low; it is also used to solve the momentum equation using the IDEAL algorithm based on the velocity pressure correction theory; The microscopic anti-diffusion model selection step is used to solve the concentration equation using a fully implicit time integration scheme; and is also used to select an anti-diffusion model to characterize the concentration distribution; The general component step includes an algebraic equation group step, a control equation discrete coefficient step, a grid coordinate step, and a coefficient zeroing step, wherein: The algebraic equation group step is used to solve the discrete momentum, temperature and concentration control equations, and is used to be called by the solid-liquid interface tracking step and the microscopic anti-diffusion model selection step; The control equation discrete coefficient step is used to solve the coefficients of the velocity and pressure equations in the momentum equation, thereby providing input for the algebraic equation system step; The grid coordinate step is used to divide the calculation area and define the grid coordinate values ​​in the calculation area; The coefficient zeroing step is used to zero the coefficients of the control equation before the control equation discrete coefficient step, so that the control equation discrete coefficient step can be repeatedly executed; Output step: used to obtain output results according to the output parameters of the brine film icing numerical model; the output results include the growth rate of the surface ice layer, ice accumulation characteristics and the influence of different factors on the ice accumulation characteristics.

[0031] It should be noted that the method described in this embodiment is applicable to the scene of seawater droplet freezing. The main solute component in seawater is sodium chloride, so the salt described in this embodiment specifically refers to sodium chloride.

[0032] In this embodiment, the typical component characteristics of polar navigation ships are structural characteristics (features) of simplified models of typical open-air equipment in ship superstructures; Typical open-air equipment includes shipboard antennas, winches, and radars; The simplified models of the typical open-air equipment include a cylindrical model, a square plate model, a thin line model, and the like.

[0033] In this implementation, the navigation environment information includes hydrological information, meteorological information and ship navigation information, wherein: Hydrological information includes seawater salinity and temperature information in the navigation area; seawater salinity information refers to the sodium chloride concentration in seawater, which is measured by a salinometer; Meteorological information includes air temperature, wind speed, humidity, droplet content and droplet diameter information in the navigation area; droplet content refers to the mass ratio of droplets per unit volume of air; Ship navigation information includes ship heading and speed information.

[0034] In this embodiment, the seawater density in the input parameter is the density of the seawater in the ship's navigation area.

[0035] In this embodiment, the seawater salinity in the input parameter is the salinity of the seawater in the ship's navigation area.

[0036] In this implementation manner, the liquid water content in the input parameter is the mass of liquid water per unit volume of air.

[0037] In this embodiment, the dynamic viscosity in the input parameters includes ice phase dynamic viscosity and liquid phase dynamic viscosity, which are used to characterize the proportionality coefficient between the friction resistance in the fluid and the relative velocity between two layers of fluid per unit distance.

[0038] In this embodiment, the initial temperature in the input parameters includes the air temperature, the seawater temperature and the surface temperature of the ship superstructure components in the initial state.

[0039] In this embodiment, the eutectic concentration in the input parameters is the minimum concentration at which the salt solution can simultaneously crystallize two solid phases with different compositions and different crystal structures at the eutectic temperature, and the eutectic concentration of the sodium chloride aqueous solution is 23.7%.

[0040] In this embodiment, the eutectic temperature in the input parameters is the lowest temperature at which the solid phases of all components in a multi-component system reach equilibrium at the same time.

[0041] In this embodiment, the partition coefficient in the input parameter is a ratio of solid-liquid concentration in a binary solidification / melting system, or can be regarded as the ratio of the slope of the liquidus line or the solidus line in a binary solidification / melting phase diagram.

[0042] In this embodiment, the liquid phase latent heat in the input parameter is the heat released when water turns into ice.

[0043] In this embodiment, the specific heat capacity in the input parameters is the amount of heat absorbed or released when a unit mass of an object changes unit temperature.

[0044] In this embodiment, the thermal diffusion coefficient in the input parameter is a physical quantity that characterizes the degree of heat diffusion in the liquid during the solidification or melting process.

[0045] In this embodiment, the component diffusion coefficient in the input parameter is a physical quantity that characterizes the degree of salt diffusion in the liquid during solidification or melting.

[0046] In this embodiment, the expansion coefficient in the input parameters includes the thermal expansion coefficient and the solute expansion coefficient, which is a physical quantity that characterizes the heat or salinity expansion capacity of the liquid during the solidification or melting process.

[0047] In this embodiment, the anti-diffusion coefficient in the input parameter is used to characterize the diffusion (anti-diffusion) degree of salt in the solid phase during solidification or melting, and its value is used to determine the type of anti-diffusion model.

[0048] In this embodiment, the mushy coefficient in the input parameter is the coefficient before the term characterizing the speed of liquid phase movement in the momentum equation in the enthalpy-porosity method, and its value range is between 104-107, and the specific value is related to the degree of supercooling.

[0049] In this embodiment, the latent heat update relaxation factor in the input parameters is used to characterize the slowdown of the update iteration from the old latent heat value to the new latent heat value in the latent heat update expression in the enthalpy-porosity method.

[0050] In this embodiment, the liquid phase fraction in the input parameters represents the solid-liquid phase state of each grid node. The liquid phase fraction value of the solid phase grid node is 0, the liquid phase fraction of the liquid phase grid node is 1, and the liquid phase fraction range of the solid-liquid mixed state grid node is (0,1).

[0051] In this implementation manner, the calculation area in the calculation parameters includes a cold plate area and an air area. The cold plate area is placed in the air area. The two areas share a set of grid systems to improve calculation efficiency.

[0052] In this implementation manner, the number of grids in the calculation parameters is the number of subdivided grids in the calculation area, which is used for control equations and main variables in the discrete numerical calculation process.

[0053] In this implementation, the time step in the calculation parameters is the time interval between each time layer in the numerical calculation process.

[0054] In this embodiment, the initial conditions in the calculation parameters are the conditions at the initial moment of the numerical calculation, including the temperature (subcooling degree), velocity and liquid fraction of the air region and the cold plate region, respectively.

[0055] In this embodiment, the boundary conditions in the calculation parameters are the conditions at the boundaries of the calculation domain (including the outer boundaries of the air area and the cold plate area), which are used to specify variables such as temperature, velocity, liquid phase fraction, contact angle, etc.; among which, the contact angle boundary conditions are applied using a dynamic contact angle model.

[0056] In this embodiment, the method involves the brine film runoff and ice growth behavior, which is a typical phenomenon of gas-liquid two-phase flow and solid-liquid phase change coupling. The studied fluid is a Newtonian fluid with the characteristics of incompressibility, laminar flow, and transient state, and the density is regarded as not changing with time, then the mass conservation equation is:

[0057] In the formula and Both represent divergence.

[0058] On this basis, the momentum equation (Navier-Stokes equation, momentum conservation equation, or NS equation for short) is obtained:

[0059] In the formula, is the density, For time, and The velocity vectors are exist and The weight in direction, is the dynamic viscosity, is the pressure term, is the gradient symbol; is the surface tension source term, which is used to reflect the surface tension effect at the gas-liquid interface; S For other source items.

[0060] In this implementation, the energy equation is expressed as follows:

[0061] In the formula, is the thermal conductivity, in units of , is the latent heat of the mixture, which can be expressed as: , is the latent heat of liquid, is the liquid phase fraction.

[0062] In this embodiment, the concentration equation is:

[0063] In the formula, is the mixed volume concentration in the solution, is the mixed mass concentration, is the diffusion coefficient of the corresponding component, is the mass source term.

[0064] In this embodiment, the temperature in the output parameter is the temperature value of each grid node in the calculation area, which is used to reflect the temperature gradient of the calculation domain to determine whether each grid node reaches the solidification temperature condition.

[0065] In this embodiment, the salinity in the output parameter is the salinity value of each grid node in the calculation area, which is used to reflect the concentration (content) of sodium chloride in the water in the ice water pool and reflect the salt precipitation during the freezing process.

[0066] In this implementation manner, the speed in the output parameter is used to characterize the speed of fluid flow.

[0067] In this implementation manner, the pressure in the output parameter is used to characterize the force acting on the contact surfaces of two objects.

[0068] In this embodiment, the liquid phase fraction in the output parameter is a basic parameter for distinguishing between solid and liquid phases in the enthalpy-porosity method, representing the solid-liquid phase state of each grid node. The liquid phase fraction value of the solid phase grid node is 0, the liquid phase fraction of the liquid phase grid node is 1, and the liquid phase fraction value range of the solid-liquid mixed state grid node is (0,1).

[0069] In this embodiment, the volume fraction in the output parameter is a basic parameter in the VOSET method, which is used to characterize the gas-liquid two-phase distribution: If the volume fraction value of a grid node is 1, the fluid in the grid node is liquid water; If the volume fraction value is 0, the fluid in the grid node is gas; If the volume fraction value is between 0 and 1, the fluid in the grid node is a gas-liquid mixture; 0.5 corresponds to the gas-liquid interface.

[0070] In this embodiment, the signed distance function in the output parameters is a basic parameter in the VOSET method, which is used to cooperate with the volume fraction to improve the solution of physical quantities at the gas-liquid interface; it is also used to directly represent the distribution of the gas-liquid interface. When its value is positive, the fluid in the corresponding grid is liquid phase, when it is negative, the fluid in the corresponding grid is gas phase, and when its value is zero, the corresponding grid is a gas-liquid interface grid.

[0071] In this embodiment, the value of the mixture fraction in the output parameter is half of the sum of the liquid fraction and the volume fraction, which is used to solve the problem that the liquid fraction or the volume fraction alone cannot simultaneously reflect and distinguish the three phases of solid, liquid and gas.

[0072] The mixture fraction is a basic parameter introduced in this implementation mode, which is an improvement on the prior art and has originality.

[0073] In this embodiment, the volume fraction is used to describe the gas-liquid interface, and the liquid phase fraction is used to describe the solid-liquid interface. Both fractions can only display the two-phase interface, and cannot simultaneously display the gas-liquid-solid interface. For this reason, a mixed fraction is introduced, and its value is equal to half of the sum of the volume fraction and the liquid phase fraction. In this way, the three phases of air-water-ice in the brine film freezing process can be displayed simultaneously.

[0074] It should be noted that in order to simultaneously reflect the gas-liquid-solid three-phase interface, the introduction of the mixture fraction comprehensively considers the influence of the two fractions. To distinguish the gas, liquid and solid three phases, as shown in Table 1.

[0075] Table 1 Three volume fraction values

[0076] In this embodiment, the surface ice layer growth rate in the output result is obtained by extracting the ice layer thickness value in the mixing fraction cloud map at different times.

[0077] In this embodiment, the ice accumulation characteristics in the output result include dense ice and porous ice, where: Dense ice has few impurities or bubbles in the ice body, and is smooth and transparent; Porous ice is ice with a large number of bubbles distributed inside the ice body, which appears milky white.

[0078] In this embodiment, the influence of different factors in the output results on the ice accumulation characteristics includes: the influence of typical component characteristics of polar navigation ships and navigation environment information on the ice accumulation characteristics.

[0079] The characteristics of typical components of polar navigation ships are the characteristics of simplified models of typical open-air equipment such as shipboard antennas, capstans, radars, etc. in the ship's superstructure, including cylinders, square plates, thin wires, etc.; Navigation environment information includes hydrological information, meteorological information and ship navigation information, including: Hydrological information includes seawater salinity and temperature in the navigation area. Salinity refers to the concentration of sodium chloride in seawater, which is measured by a salinometer. Meteorological information includes air temperature, wind speed, humidity, droplet content, droplet diameter, etc. in the navigation area, where droplet content refers to the mass ratio of droplets per unit volume of air; The ship navigation information includes the ship's heading, speed, etc.

[0080] In addition, in one embodiment, in the grid coordinate step: A structured staggered grid system is established using the finite volume method to make the velocity solution more stable.

[0081] In this implementation, a structured staggered grid system is established through a finite volume method to make the velocity solution more stable.

[0082] In this embodiment, the finite volume method and staggered grids are used to discretize the basic control equations.

[0083] In this embodiment, a structured staggered grid is adopted, and three sets of grid schemes are used in the two-dimensional grid, which can effectively avoid the divergence of pressure gradient.

[0084] In addition, in one embodiment, the cold plate surface contact angle boundary condition is set to reflect the different wettability of the cold plate surface and consider the surface hysteresis effect caused by the different wettability, including: The dynamic contact angle of the water droplet during its motion is solved in real time, and the dynamic contact angle input is provided for the gas-liquid interface tracking step, thereby correcting the signed distance function value at the solid boundary, thereby reflecting the different wettability of the cold plate surface in the numerical model of brine film icing; The different wettabilities include hydrophilicity, hydrophobicity and superhydrophobicity; The hydrophilicity is a contact angle of less than 90°; The hydrophobicity is a contact angle between 90° and 150°; The super hydrophobicity is a contact angle greater than 150°; The surface hysteresis effect caused by different wettabilities is the dynamic change of the contact angle between the water droplet and the solid surface during the spreading and retracting process.

[0085] It should be noted that the spreading of water droplets on solid surfaces is a typical wetting phenomenon. Its ability to spread on solid surfaces is called wettability, which is affected by the roughness of the solid surface (roughness refers to the unevenness composed of small spacing and tiny peaks and valleys on the processed surface). Wettability affects the spreading and shrinking process of water droplets on the surface, and its strength is usually quantitatively expressed by the contact angle. During the movement of water droplets on the wall, the contact line of the gas-liquid interface at the air-liquid-wall three-phase contact point will form a certain angle with the wall. The angle close to the liquid phase is called the contact angle. The contact angle has an important influence on the dynamic behaviors of water droplets such as spreading, retraction, and deformation. According to the size of the contact angle, wettability can be divided into several types, such as hydrophilic (contact angle less than 90°), hydrophobic (contact angle between 90° and 150°), and superhydrophobic (contact angle greater than 150°).

[0086] When a water drop remains stationary on a smooth surface, the corresponding contact angle is usually called the static contact angle. .

[0087] During the movement of water droplets, on the one hand, they move forward due to the initial kinetic energy, and on the other hand, they are hindered by the hysteresis caused by the surface roughness. Due to the opposing effects of kinetic energy and hysteresis effect, the contact angle between the water droplet and the surface has the characteristic of dynamic change. This changing contact angle is called dynamic contact angle. .

[0088] The dynamic contact angle is usually around the static contact angle Fluctuation, greater than Maximum contact angle Advancing contact angle , less than The minimum contact angle is the receding contact angle .

[0089] In addition, in one embodiment, the real-time solution of the dynamic contact angle during the motion of the water droplet includes: The Kistler dynamic contact angle model is used to solve the dynamic contact angle. :

[0090] In the formula, is the Hoffman smoothing function, which is calculated as follows:

[0091] In the formula, the independent variable Defined as:

[0092] In the formula, is the inverse function of the Hoffman smoothing function, which is solved by the following iteration:

[0093] in, is the capillary number, which is used to characterize the ratio of viscosity to surface tension; To balance the contact angle, The inverse function of , used to reflect the contact angle hysteresis effect; According to the three-dimensional three-phase contact point velocity of water droplet-air-fixed surface direction to determine :

[0094] In the actual scene of water drop movement, the following formula is used to solve :

[0095] in,

[0096] in, is the distance between the three-phase contact point and the initial impact point of the water droplet.

[0097] It should be noted that the capillary number Three-phase contact speed About, when Sometimes, there are , which is equivalent to a fixed dynamic contact angle value that does not change; however, in the actual scene of water droplet movement, the capillary number will not be always 0, so the following formula is used in this implementation to solve : .

[0098] It should be noted that In order to balance the contact angle, a fixed value is set, which is equal to the advancing contact angle or the receding contact angle. This value is introduced to reflect the contact angle hysteresis effect (i.e., surface hysteresis effect).

[0099] It should be noted that It is a dynamic contact angle, changing in real time.

[0100] It should be noted that It is the static contact angle, which is fixed.

[0101] It should be noted that the dynamic contact angle Capillary number and equilibrium contact angle about, and and The three-phase contact speed Impact, that is The accuracy of the solution has a great influence on the dynamic contact angle.

[0102] Additionally, in one embodiment: The solution is obtained by the following method: By calculating the coordinates of the grid nodes where the three-phase contact points are located (which can be converted into ) to obtain :

[0103] In the formula, is the distance between the three-phase contact point and the initial impact point of the water drop; The above formula is generally discretized into the following form for solution:

[0104] In the formula, is the lateral distance between the three-phase contact point of the current layer water droplet and the initial impact point, represent The horizontal distance value before time; Equal to 1 or more time steps .

[0105] visible, The selection of has a great influence on the calculation results.

[0106] In this embodiment, It is obtained by transforming the coordinate values ​​of the grid nodes where the three-phase contact points are located.

[0107] It should be noted that the contact line solution method is improved in this embodiment, so that the obtained real-time contact angle is more continuous and accurate.

[0108] In addition, in one embodiment: according to the change of the spreading front, the Size value: If the current spreading front position does not change (i.e. unchanged), then Take the initial value or the value of the previous period; If the spreading front position changes, then , the corresponding distance mark of the current three-phase contact point transverse coordinate is ; Calculate according to the following formula :

[0109] In the formula, is the time corresponding to the current time layer, for Corresponding time; the superscript new indicates the current time layer, and old indicates the value before Dt time; Calculation completed After that, , ; Get accurate three-phase contact point speed After that, the dynamic contact angle can be obtained based on the Kistler dynamic contact angle model. .

[0110] It should be noted that Figure 2 It shows a schematic diagram of the implementation process of the dynamic contact angle application step (component); this step (component) is preceded by an "input step (module)" to obtain input parameters; after the dynamic contact angle is calculated, an influence is exerted on the gas-liquid interface tracking step (component) and the isolated cold plate step (component).

[0111] It should be noted that the air-liquid interface tracking step (component) is used to determine the air-liquid interface during the motion of the water droplet, and the parameters to be solved include the signed distance function and volume fraction. Among them, the behavior of water droplets spreading and retracting on the surface is affected by the signed distance function value at the boundary At the same time, the dynamic contact angle exerts an influence on the signed distance function at the boundary, thereby affecting the motion of the water droplet, thus reflecting the influence of the dynamic contact angle on the motion of the water droplet. Therefore, the influence of the contact angle is reflected in the signed distance function value at the boundary of the solid surface. On the solution.

[0112] In addition, in one embodiment: the signed distance function within the grid at the boundary The solution process is as follows: Take the case of water droplets spreading to the right on a horizontal plane as an example: The grid node where the three-phase contact point is located is ; and Node and its adjacent inner nodes The distance from the gas-liquid interface, then the angle between the gas-liquid interface and the horizontal wall The following conditions must be met:

[0113] When the liquid phase is the main phase, the liquid region If is a negative value, the above formula is transformed into programming form:

[0114] From the above formula, we can see that the signed distance function in the solid boundary grid is The solution of needs to determine the dynamic contact angle first; When executing the gas-liquid interface tracking step to track the gas-liquid interface, obtain the nodes in the current time layer calculation area. back, Call the dynamic contact angle application step (component) to obtain the dynamic contact angle value, and then calculate the outermost node through the above formula In this way, the influence of the dynamic contact angle is successfully exerted by modifying the value of the signed distance function within the solid boundary mesh.

[0115] Signed distance function within the grid at the boundary The solution process, such as Figure 3 shown.

[0116] In addition, in one embodiment, in the isolated cold plate step, the isolated cold plate (solid) is incorporated into the fluid calculation domain and solved as a whole; the overall processing idea is to make the velocity in the isolated cold plate area zero and the temperature a given sub-zero value. : The large coefficient method is used to make the temperature in the isolated cold plate area take a given value , that is, before solving the temperature algebraic equations, let the nodes in the isolated cold plate region The corresponding discrete equation is: ; In the formula, , ,in, is a larger value (such as ); This is done to isolate the nodes in the cold plate area. Corresponding temperature Always take the given value ; The condition that the velocity field in the isolated cold plate area is always zero can be achieved by imitating the above method: The corresponding position value is 0; However, the velocity value will be corrected during the pressure-velocity coupling solution process. To avoid the pressure correction of the velocity in the isolated cold plate area, thereby destroying the zero velocity condition, the following settings need to be repeated in each time step: The first setting is: before each iteration, the velocity in the isolated cold plate area is set to 0. Since the staggered grid is used, when For isolated cold plate area boundary nodes, additional , similarly When is the right boundary, ; The second setting: Before solving the velocity algebraic equations, let the main diagonal coefficients of each velocity discrete matrix in the isolated cold plate region be is a larger value (such as ), thereby ensuring the calculated intermediate speed ; Setting 3: Use the IDEAL algorithm to solve the pressure-velocity coupling. The pressure correction value needs to be calculated in two inner iterations to avoid the pressure correction value in the isolated cold plate area affecting the velocity. To make corrections, we need to set the speed correction equation coefficient ( etc.) to take a minimum value (such as ).

[0117] In this embodiment, an isolated cold plate numerical method is used to achieve unified solution of the fluid domain and the solid domain in the same grid system, avoiding separate modeling of the fluid grid and the solid grid, which is easier to implement through programming.

[0118] In addition, in one embodiment, in the solid-liquid interface tracking step, the temperature and liquid fraction of each grid node in the calculation area are obtained by the following method: The source term of the energy equation contains a transient term of the liquid fraction (i.e., in the process of solving the energy equation, the current time layer contains two unknown quantities: temperature and liquid fraction); Iteratively solve the liquid fraction and obtain the temperature after the convergence result is obtained; When solving iteratively, only one iteration is required at a time: After the solution of the layer velocity and pressure at the current time is completed, the liquid phase fraction is updated and its upper and lower limits are specified, thereby obtaining the source term of the energy equation and the coefficient of the discrete equation; Using the steps of algebraic equations, solve the energy equation to obtain the temperature value, determine the source term and the discrete equation coefficient, and solve the algebraic equations to obtain the current temperature; The above steps are iterated several times until the temperature and liquid fraction of the current time layer reach consistent convergence conditions.

[0119] In this implementation, both the convection term and the diffusion term of the control discrete equation are in a power format.

[0120] It should be noted that the source term of the energy equation contains the transient term of the liquid fraction, which means that the source term of the energy equation contains both the liquid fraction value at the previous moment and the liquid fraction value at the current moment, that is, there are two unknown quantities (temperature and liquid fraction) in the energy equation of the current time layer. If the liquid fraction is updated directly using the piecewise relationship between the liquid fraction and the temperature, it is easy to cause calculation divergence.

[0121] In this implementation, in order to improve convergence and avoid computational divergence, the liquid fraction is iteratively solved to obtain a liquid fraction value consistent with the temperature solution at the current moment.

[0122] It should be noted that the existing liquid fraction iteration method performs two iterations in the same time layer, which has low computational efficiency.

[0123] In this implementation, only one iteration is required in the same layer, and the calculation efficiency is high.

[0124] In addition, in one embodiment, the volume fraction is used to characterize the gas-liquid two-phase distribution: If the volume fraction value of a grid node is 1, the fluid in the grid node is liquid water; If the volume fraction value is 0, the fluid in the grid node is gas; If the volume fraction value is between 0 and 1, the fluid in the grid node is a gas-liquid mixture, where 0.5 corresponds to the gas-liquid interface.

[0125] In this embodiment, the signed distance function is used to directly represent the distribution of the gas-liquid interface: When the value of the signed distance function is positive, the fluid in the corresponding grid is in liquid phase; When the value of the signed distance function is negative, the fluid in the corresponding grid is gas phase. When the value of the signed distance function is zero, the corresponding grid is a gas-liquid interface grid.

[0126] In addition, in one embodiment, the coefficient of the source term is set to maintain fluidity when the liquid phase ratio is high and to suppress flow when the liquid phase ratio is low: When the droplet region is entirely liquid, the source term is 0, and the original two-phase flow momentum equation is not affected; When the droplet begins to freeze, the interior of the droplet becomes an ice-water mixed region. At this time, the liquid phase fraction is between (0,1) and gradually decreases. As a result, the source term accounts for an increasingly larger proportion in the momentum equation and gradually exceeds the transient term, convection term, and diffusion term. Until the droplet is completely transformed into the solid phase, the source term reaches a maximum value, so the velocity of the solid phase grid node is basically 0, realizing the transition from liquid phase to solid phase.

[0127] It should be noted that the momentum equation used in the traditional VOSET method only contains commonly used transient terms, convection terms, diffusion terms, surface tension terms, and body terms. In this embodiment, while considering gas-liquid two-phase flow phenomena such as air flow and liquid movement, the freezing of the free liquid surface in the ice-water object is also considered. For this reason, the source term about the liquid phase fraction is added to the momentum equation in this embodiment to reflect the effect of freezing on movement.

[0128] It should be noted that one of the technical difficulties involved in the simulation of the brine film freezing process is the coupled tracking of the gas-liquid two-phase flow interface such as cold air flow and water surface fluctuations and the solid-liquid two-phase phase change interface such as water surface freezing. The gas-liquid two-phase flow interface focuses on solving the momentum equation and updating the fluid volume fraction, while the water surface freezing interface focuses on solving the energy equation and updating the liquid phase fraction. How to comprehensively consider and solve the speed, pressure, temperature and other parameters for solving the two types of equations is a major difficulty.

[0129] In this implementation, the momentum equation is solved mainly by using the IDEAL algorithm based on the velocity-pressure correction theory, which can improve the speed of solving the velocity-pressure.

[0130] In this embodiment, based on the IDEAL algorithm, an energy equation is added to solve, and a source term about the liquid phase fraction is added to the energy equation, and a source term about the liquid phase fraction is also added to the momentum equation, thereby establishing a connection between the energy equation and the momentum equation.

[0131] In this embodiment, in order to reflect the effect of freezing on liquid flow, the relationship between the momentum equation and the energy equation is improved, and the liquid fraction Directly reflects the solid-liquid mixing state in the mushy region and affects the two-phase flow during the phase change process. This effect can be achieved by adding source terms to the basic NS equations. (Assuming that the motion in the mushy region obeys Darcy's law for flow in porous media and conforms to the Carman-Koseny hypothesis), It can be understood as the flow resistance caused by solidification, and its specific form is: ; In the formula It is a minimum value. To ensure that the denominator is not 0, it is taken as 0.001; It is a constant value in the pasty region, which can control the speed reduction rate. It is recommended to be between 104 and 107. The appropriate value can maintain fluidity when the liquid phase ratio is high and inhibit flow when the liquid phase ratio is low.

[0132] When the droplet area is entirely liquid, , i.e. the source term , the original two-phase flow NS equation is not affected; when the droplet begins to freeze, the inside of the droplet becomes an ice-water mixed area, that is, a mushy area. Between (0,1) and gradually decreases, thus It accounts for an increasingly larger proportion in the NS equations and gradually exceeds the transient term, convection term and diffusion term until the droplet is completely transformed into a solid phase and the velocity is basically zero.

[0133] It should be noted that another technical difficulty involved in the simulation of the brine membrane freezing process is how to consider the influence of salinity. Under the combined action of concentration (salinity) and temperature gradient, solutes will precipitate and move toward the liquid phase in the fluid in the ice-water mixing area, forming brine channels in the process. That is, brine channels are mainly caused by heat-salt double diffusion convection. Compared with the natural convection caused by temperature gradient in freshwater freezing, the heat-solute double diffusion convection in brine freezing is more complicated. Its mechanism of action on the formation and evolution of brine channels is still unclear, and no relevant literature has been found. It is still a major problem that needs to be solved in the study of seawater freezing. The key to solving this problem is to accurately track the various interfaces in the heat-solute double diffusion convection process in the brine membrane, and reflect the coupling effect of interface dynamics and heat-solute driving force.

[0134] It should be noted that, in order to solve the problem that the existing simulation of the brine membrane freezing process cannot take into account the dynamic influence of salinity, this embodiment introduces a concentration equation while solving the momentum equation and the energy equation.

[0135] Considering that the main solute component in seawater is sodium chloride, this embodiment regards seawater as a binary mixture of sodium chloride and pure water. On the basis of the pure water freezing model, the salinity (concentration) equation is introduced. Based on the local thermodynamic equilibrium assumption and phase diagram, the solidus and liquidus temperatures in the energy equation are determined by the salinity, thereby realizing the coupling of the temperature field and the solute field.

[0136] At the same time, the body and source terms of the momentum equation are improved to reflect the effects of thermal buoyancy, solute buoyancy and forced convection on salt diffusion and convection, and to simulate the heat-solute double diffusion convection process under the combined action of temperature and salinity gradients.

[0137] The freezing of brine film includes not only macroscopic phenomena such as heat-solute double diffusion convection and ice-water interface advancement, but also local segregation phenomena such as salt expulsion from dendrites and back diffusion. Therefore, it is necessary to couple the analysis of macroscopic and local segregation processes. It is proposed to establish a back diffusion model based on the lever principle and the Scheil-Gulliver principle, correct the source term of the macroscopic salinity equation, describe the role of salinity in controlling phase change, and reflect the redistribution of salt through diffusion transport.

[0138] It should be noted that traditional velocity and pressure correction methods include SIMPLE, SIMPLER, SIMPLEC, PISO and other algorithms, most of which require assumptions. In view of the gas-liquid-solid three-phase coupling tracking problem, this implementation first uses the pressure correction IDEAL algorithm to improve the convergence efficiency of velocity and pressure solutions, thereby constructing a fluid calculation domain, which can overcome the two basic assumptions in the traditional method, making the calculation process more efficient and stable, and on this basis, the VOSET method is used to achieve gas-liquid interface tracking.

[0139] It should be noted that the mainstream methods for tracking the gas-liquid interface are VOF and Level Set methods. The former has the advantage of mass conservation, but the disadvantage is that the physical quantities at the interface are discontinuous, while the latter has the advantage of smooth and continuous physical quantities at the interface, but the disadvantage is that it cannot meet the mass conservation. This implementation adopts the VOSET method to track the gas-liquid interface, combining the advantages of the VOF method and the Level Set method, which can not only conserve mass but also ensure the smooth and continuous parameters such as density and surface tension at the interface.

[0140] It should be noted that although the VOSET method and the enthalpy-porosity method are used for tracking gas-liquid interface and solid-liquid interface respectively in existing methods (such as water droplet impact spreading dynamics, metal melting, welding, etc.), there is no literature combining the two to track gas-liquid-solid three-phase, and there is no literature combining the two to simulate the freezing process of brine membranes. The difficulty lies in that the two methods have different focuses, and it is difficult to establish the relationship between the parameters, such as the effect of temperature on the liquid fraction. , while the velocity affects the volume fraction change , it is difficult to explain how the two scores affect each other, and only one score can be displayed in the same cloud chart.

[0141] In this implementation, a fully implicit time integration scheme is used to solve the concentration equation to ensure unconditional stability of the result.

[0142] Liquid concentration in a single grid Uniform distribution, while the diffusion of solute in the solid phase is relatively slow, which will form a solute gradient, thus causing solute back diffusion. The solid phase concentration in a single grid It is not evenly distributed and requires points to obtain.

[0143] Brine film freezing is a non-equilibrium solidification, and a reasonable reverse diffusion model needs to be selected to determine the liquid phase concentration. , thereby correcting the latent heat update method in the source term of the energy equation in the pure water solidification model (i.e., selecting the anti-diffusion model to characterize the concentration distribution); among them, the anti-diffusion model includes the lever principle, the Scheil-Gulliver principle, the Voller improved anti-diffusion model and other models.

[0144] Additionally, in one embodiment, selecting a reverse diffusion model to characterize the concentration distribution includes: In the pure water solidification model, the inverse function in the latent heat update expression is is a fixed value; After considering salt segregation, the latent heat update expression is improved to:

[0145] In the formula, It is not a constant value, but a value that changes with the mixture concentration during the phase change process; thus, the dynamic relationship between the temperature field and the solute field is established; After further considering the solute back diffusion, the latent heat update expression is improved to a specific form based on each back diffusion model, and the Scheil-Gulliver non-equilibrium theory is substituted into it, then:

[0146] In the Leverage Principle:

[0147] In Voller's improved anti-diffusion model, which is between Scheil-Gulliver non-equilibrium theory and the lever principle: .

[0148] The above-mentioned reverse diffusion model is selected to characterize the concentration distribution and a numerical model of brine film freezing is constructed.

[0149] It should be noted that the numerical model of brine film freezing includes a dynamic contact angle application step (component), an isolated cold plate step (component), a solid-liquid interface tracking step (component), a gas-liquid interface tracking step (component), a microscopic anti-diffusion model selection step (component) and a general step (component). Among them, the microscopic anti-diffusion model selection step (component) is used to select the anti-diffusion model to characterize the concentration distribution, which reflects the salt precipitation process during the brine film freezing process.

[0150] In addition, in one embodiment, a device for simulating the freezing process of salt water film on the surface of typical components of polar ships is provided: The device is used to implement the above-mentioned method for simulating the freezing process of salt water film on the surface of typical components of polar ships.

[0151] The device comprises the following modules: Input module: used to obtain the input parameters of the salt water film icing numerical model according to the characteristics of typical components of polar ships and navigation environment information; The input parameters of the numerical model of brine film icing include physical parameters and calculation parameters; The physical property parameters include seawater density and salinity, liquid water content, dynamic viscosity, initial temperature, eutectic concentration, eutectic temperature, partition coefficient, liquid phase latent heat, specific heat capacity, thermal diffusion coefficient, component diffusion coefficient, expansion coefficient, back diffusion coefficient, pastiness coefficient, latent heat renewal relaxation factor and liquid phase fraction; The calculation parameters include calculation area, number of grids, time step, initial conditions, boundary conditions, water drop release interval and number of water drops; Solver: used to solve the momentum equation, energy equation and concentration equation according to the input parameters of the brine film icing numerical model, and obtain the output parameters of the brine film icing numerical model; The momentum equation is a fluid equation established in accordance with Newton's second law; The energy equation is an energy conservation equation, which complies with the first law of thermodynamics and takes temperature as the variable to be determined; The concentration equation is a component conservation equation that complies with the component mass conservation law and is used to obtain the concentration parameter value in the calculation domain; The output parameters are temperature, salinity, velocity, pressure, liquid fraction, volume fraction, signed distance function and mixing fraction; The solution module includes a dynamic contact angle application component, an isolated cold plate component, a solid-liquid interface tracking component, a gas-liquid interface tracking component, a microscopic anti-diffusion model selection component and a general component; The dynamic contact angle applying component is used to set the contact angle boundary condition of the cold plate surface to reflect the different wettability of the cold plate surface and consider the surface hysteresis effect caused by the different wettability; The isolated cold plate component is used to set an isolated cold plate area that is not connected to the outer boundary in the air calculation domain to simulate the runoff phenomenon of the salt water film; The solid-liquid interface tracking component is used to solve the energy equation and obtain the temperature and liquid fraction of each grid node in the calculation area; The gas-liquid interface tracking component is used to solve the volume fraction and signed distance function that characterize the gas-liquid interface, and selects the VOSET gas-liquid interface tracking method to process the gas-liquid interface changes during the motion of the water droplet; it is also used to add a source term about the liquid phase fraction to the momentum equation, and by setting the coefficient of the source term, the fluidity is maintained when the liquid phase ratio is high, and the flow is suppressed when the liquid phase ratio is low; it is also used to solve the momentum equation using the IDEAL algorithm based on the velocity pressure correction theory; The microscopic anti-diffusion model selection component is used to solve the concentration equation using a fully implicit time integration scheme; it is also used to select an anti-diffusion model to characterize the concentration distribution; The general component includes an algebraic equation group module, a control equation discrete coefficient module, a grid coordinate module and a coefficient zeroing module, wherein: The algebraic equation group module is used to solve the discrete momentum, temperature and concentration control equations, and is used to be called by the solid-liquid interface tracking step and the microscopic anti-diffusion model selection step; The control equation discrete coefficient module is used to solve the coefficients of the velocity and pressure equations in the momentum equation, thereby providing input for the algebraic equation system step; The grid coordinate module is used to divide the calculation area and define the grid coordinate values ​​in the calculation area; The coefficient zeroing module is used to zero the coefficients of the control equation before the discrete coefficient step of the control equation, so that the discrete coefficient step of the control equation can be repeatedly executed; Output module: used to obtain output results according to the output parameters of the brine film icing numerical model; the output results include the growth rate of the surface ice layer, ice accumulation characteristics and the influence of different factors on the ice accumulation characteristics.

[0152] The technical solution provided by the present invention is further described in detail above through several specific implementation modes in order to highlight the advantages and benefits of the technical solution provided by the present invention. However, the several specific implementation modes described above are not intended to be used as limitations on the present invention. Any reasonable changes and improvements to the present invention, reasonable combinations of implementation modes and equivalent substitutions based on the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for simulating the freezing process of salt water film on the surface of typical components of polar ships, characterized in that: The method comprises the following steps: Input step: used to obtain input parameters of the salt water film icing numerical model according to the characteristics of typical components of polar navigation ships and navigation environment information; The input parameters of the numerical model of brine film icing include physical parameters and calculation parameters; The physical property parameters include seawater density and salinity, liquid water content, dynamic viscosity, initial temperature, eutectic concentration, eutectic temperature, partition coefficient, liquid phase latent heat, specific heat capacity, thermal diffusion coefficient, component diffusion coefficient, expansion coefficient, back diffusion coefficient, pastiness coefficient, latent heat renewal relaxation factor and liquid phase fraction; The calculation parameters include calculation area, number of grids, time step, initial conditions, boundary conditions, water drop release interval and number of water drops; Solving step: solving the momentum equation, energy equation and concentration equation according to the input parameters of the brine film icing numerical model, and obtaining the output parameters of the brine film icing numerical model; The momentum equation is a fluid equation established in accordance with Newton's second law; The energy equation is an energy conservation equation, which complies with the first law of thermodynamics and takes temperature as the variable to be determined; The concentration equation is a component conservation equation that complies with the component mass conservation law and is used to obtain the concentration parameter value in the calculation domain; The output parameters are temperature, salinity, velocity, pressure, liquid fraction, volume fraction, signed distance function and mixing fraction; The solving steps include a dynamic contact angle application step, an isolated cold plate step, a solid-liquid interface tracking step, a gas-liquid interface tracking step, a microscopic anti-diffusion model selection step, and a general component step; The dynamic contact angle application step is used to set the contact angle boundary conditions on the cold plate surface to reflect the different wettabilities on the cold plate surface and to consider the surface hysteresis effect caused by the different wettabilities; The isolated cold plate step is used to set an isolated cold plate area that is not connected to the outer boundary in the air calculation domain to simulate the runoff phenomenon of the salt water film; The solid-liquid interface tracking step is used to solve the energy equation and obtain the temperature and liquid fraction of each grid node in the calculation area; The gas-liquid interface tracking step is used to solve the volume fraction and signed distance function that characterize the gas-liquid interface, and select the VOSET gas-liquid interface tracking method to process the gas-liquid interface changes during the motion of the water droplet; it is also used to add a source term about the liquid phase fraction to the momentum equation, and by setting the coefficient of the source term, the fluidity is maintained when the liquid phase ratio is high, and the flow is suppressed when the liquid phase ratio is low; it is also used to solve the momentum equation using the IDEAL algorithm based on the velocity pressure correction theory; The microscopic anti-diffusion model selection step is used to solve the concentration equation using a fully implicit time integration scheme; and is also used to select an anti-diffusion model to characterize the concentration distribution; The general component step includes an algebraic equation group step, a control equation discrete coefficient step, a grid coordinate step, and a coefficient zeroing step, wherein: The algebraic equation group step is used to solve the discrete momentum, temperature and concentration control equations, and is used to be called by the solid-liquid interface tracking step and the microscopic anti-diffusion model selection step; The control equation discrete coefficient step is used to solve the coefficients of the velocity and pressure equations in the momentum equation, thereby providing input for the algebraic equation system step; The grid coordinate step is used to divide the calculation area and define the grid coordinate values ​​in the calculation area; The coefficient zeroing step is used to zero the coefficients of the control equation before the control equation discrete coefficient step, so that the control equation discrete coefficient step can be repeatedly executed; Output step: used to obtain output results according to the output parameters of the brine film icing numerical model; the output results include the growth rate of the surface ice layer, ice accumulation characteristics and the influence of different factors on the ice accumulation characteristics.

2. The method for simulating the freezing process of salt water film on the surface of typical components of polar ships according to claim 1 is characterized in that: The value of the mixed fraction in the output parameter is half of the sum of the liquid fraction and the volume fraction, and is used to solve the problem that the liquid fraction or the volume fraction alone cannot simultaneously reflect and distinguish the three phases of solid, liquid and gas.

3. The method for simulating the freezing process of salt water film on the surface of typical components of polar ships according to claim 1 is characterized in that: The contact angle boundary condition of the cold plate surface is set to reflect the different wettability of the cold plate surface and consider the surface hysteresis effect caused by the different wettability, including: The dynamic contact angle of the water droplet during its motion is solved in real time, and the dynamic contact angle input is provided for the gas-liquid interface tracking step, thereby correcting the signed distance function value at the solid boundary, thereby reflecting the different wettability of the cold plate surface in the numerical model of brine film icing; The different wettabilities include hydrophilicity, hydrophobicity and superhydrophobicity; The hydrophilicity is a contact angle of less than 90°; The hydrophobicity is a contact angle between 90° and 150°; The super hydrophobicity is a contact angle greater than 150°; The surface hysteresis effect caused by different wettabilities is the dynamic change of the contact angle between the water droplet and the solid surface during the spreading and retracting process.

4. The method for simulating the freezing process of salt water film on the surface of typical components of polar ships according to claim 3 is characterized in that: The real-time solution of the dynamic contact angle during the water droplet movement process includes: The Kistler dynamic contact angle model is used to solve the dynamic contact angle. : In the formula, is the Hoffman smoothing function, which is calculated as follows: In the formula, the independent variable Defined as: In the formula, is the inverse function of the Hoffman smoothing function, which is solved by the following iteration: in, is the capillary number, which is used to characterize the ratio of viscosity to surface tension; To balance the contact angle, The inverse function of , used to reflect the contact angle hysteresis effect; According to the three-dimensional three-phase contact point velocity of water droplet-air-fixed surface direction to determine : In the actual scene of water drop movement, the following formula is used to solve : in, in, is the distance between the three-phase contact point and the initial impact point of the water droplet.

5. The method for simulating the freezing process of salt water film on the surface of typical components of polar ships according to claim 1, characterized in that: In the isolating cold plate step: The isolated cold plate is included in the fluid calculation domain and solved as a whole, so that the velocity in the isolated cold plate area is zero and the temperature is a given sub-zero value. : The large coefficient method is used to make the temperature in the isolated cold plate area take a given value : Before solving the temperature algebraic equations, let the nodes in the isolated cold plate region The corresponding discrete equation is: ； In the formula, , ,in, is not less than a larger value of ; To isolate the nodes in the cold plate area Corresponding temperature Always take the given value ; Will The corresponding position is set to 0, so that the velocity field in the isolated cold plate area is always zero; To avoid pressure correction to velocity in isolated cold plate regions, which would destroy the zero velocity condition, the following settings are repeated in each time step: The first setting is: before each iteration, the velocity in the isolated cold plate area is set to 0; because of the staggered grid, when When is the boundary node in the isolated cold plate area, let , similarly When is the right boundary, ; The second setting: Before solving the velocity algebraic equations, let the main diagonal coefficients of each velocity discrete matrix in the isolated cold plate region be is not less than A larger value of , thereby ensuring that the calculated intermediate speed ; The third setting: Use the IDEAL algorithm to solve the pressure-velocity coupling. Calculate the pressure correction value in two inner iterations. To avoid the pressure correction value in the isolated cold plate area to correct the velocity, set the velocity correction equation coefficient to be not greater than The minimum value of .

6. The method for simulating the freezing process of salt water film on the surface of typical components of polar ships according to claim 1, characterized in that: In the solid-liquid interface tracking step, the temperature and liquid fraction of each grid node in the calculation area are obtained by the following method: The source term of the energy equation contains a transient term of liquid fraction; Iteratively solve the liquid fraction and obtain the temperature after the convergence result is obtained; When solving iteratively, only one iteration is required at a time: After the solution of the layer velocity and pressure at the current time is completed, the liquid phase fraction is updated and its upper and lower limits are specified, thereby obtaining the source term of the energy equation and the coefficient of the discrete equation; Using the steps of algebraic equations, solve the energy equation to obtain the temperature value, determine the source term and the discrete equation coefficient, and solve the algebraic equations to obtain the current temperature; The above steps are iterated several times until the temperature and liquid fraction of the current time layer reach consistent convergence conditions.

7. A device for simulating the freezing process of salt water film on the surface of typical components of polar ships, characterized in that: The device comprises the following modules: Input module: used to obtain the input parameters of the salt water film icing numerical model according to the characteristics of typical components of polar ships and navigation environment information; The input parameters of the numerical model of brine film icing include physical parameters and calculation parameters; The physical property parameters include seawater density and salinity, liquid water content, dynamic viscosity, initial temperature, eutectic concentration, eutectic temperature, partition coefficient, liquid phase latent heat, specific heat capacity, thermal diffusion coefficient, component diffusion coefficient, expansion coefficient, back diffusion coefficient, pastiness coefficient, latent heat renewal relaxation factor and liquid phase fraction; The calculation parameters include calculation area, number of grids, time step, initial conditions, boundary conditions, water drop release interval and number of water drops; Solver: used to solve the momentum equation, energy equation and concentration equation according to the input parameters of the brine film icing numerical model, and obtain the output parameters of the brine film icing numerical model; The momentum equation is a fluid equation established in accordance with Newton's second law; The energy equation is an energy conservation equation, which complies with the first law of thermodynamics and takes temperature as the variable to be determined; The concentration equation is a component conservation equation that complies with the component mass conservation law and is used to obtain the concentration parameter value in the calculation domain; The output parameters are temperature, salinity, velocity, pressure, liquid fraction, volume fraction, signed distance function and mixing fraction; The solution module includes a dynamic contact angle application component, an isolated cold plate component, a solid-liquid interface tracking component, a gas-liquid interface tracking component, a microscopic anti-diffusion model selection component and a general component; The dynamic contact angle applying component is used to set the contact angle boundary condition of the cold plate surface to reflect the different wettability of the cold plate surface and consider the surface hysteresis effect caused by the different wettability; The isolated cold plate assembly is used to set an isolated cold plate area that is not connected to the outer boundary in the air calculation domain to simulate the runoff phenomenon of the salt water film; The solid-liquid interface tracking component is used to solve the energy equation and obtain the temperature and liquid fraction of each grid node in the calculation area; The gas-liquid interface tracking component is used to solve the volume fraction and signed distance function that characterize the gas-liquid interface, and selects the VOSET gas-liquid interface tracking method to process the gas-liquid interface changes during the motion of the water droplet; it is also used to add a source term about the liquid phase fraction to the momentum equation, and by setting the coefficient of the source term, the fluidity is maintained when the liquid phase ratio is high, and the flow is suppressed when the liquid phase ratio is low; It is also used to solve the momentum equation using the IDEAL algorithm based on velocity pressure correction theory; The microscopic anti-diffusion model selection component is used to solve the concentration equation using a fully implicit time integration scheme; it is also used to select an anti-diffusion model to characterize the concentration distribution; The general component includes an algebraic equation group module, a control equation discrete coefficient module, a grid coordinate module and a coefficient zeroing module, wherein: The algebraic equation group module is used to solve the discrete momentum, temperature and concentration control equations, and is used to be called by the solid-liquid interface tracking step and the microscopic anti-diffusion model selection step; The control equation discrete coefficient module is used to solve the coefficients of the velocity and pressure equations in the momentum equation, thereby providing input for the algebraic equation system step; The grid coordinate module is used to divide the calculation area and define the grid coordinate values ​​in the calculation area; The coefficient zeroing module is used to zero the coefficients of the control equation before the discrete coefficient step of the control equation, so that the discrete coefficient step of the control equation can be repeatedly executed; Output module: used to obtain output results according to the output parameters of the brine film icing numerical model; the output results include the growth rate of the surface ice layer, ice accumulation characteristics and the influence of different factors on the ice accumulation characteristics.

8. A computer device comprising: A processor and a memory, characterized in that the memory is used to store executable instructions of the processor, and the processor is configured to execute the method for simulating the icing process of salt water film on the surface of typical components of polar ships as described in any one of claims 1 to 6 by executing the executable instructions.

9. A computer storage medium, characterized in that: The storage medium stores a computer program, and when the computer program is run, the method for simulating the freezing process of salt water film on the surface of typical components of polar ships according to any one of claims 1 to 6 is executed.

10. A computer program product comprising a computer program / instructions, characterized in that When the computer program / instructions are executed by a processor, the steps of the method for simulating the freezing process of salt water film on the surface of typical components of polar ships as described in any one of claims 1 to 6 are implemented.

Images