Towards the efficient modelling of trapped air pockets during squeeze flow

In most bonding processes, an adhesive is applied to a substrate in a specific pattern before the second substrate is subsequently pressed against it. During this, the adhesive flows in such a way that, ideally, it completely fills the joint. In practice, however, areas with entrapped air frequently remain in the bonded adhesive layer. Within the scope of a research project, these flows are systematically analyzed in order to identify optimal initial application patterns for the adhesive and substrate geometry to minimise such risks. For this purpose, the authors use an efficient flow model, the partially filled gaps model (PFGM), extended in this study to include the functionality of trapped air pockets. Depending on the volume fractions of air and adhesive, the flow of both phases is computed. Therefore, the model is introduced and fully described, benchmarked with respect to its plausibility and functionality, and results obtained are compared with a CFD calculation. Thereafter, the functionality of openings and closings of the pockets are analyzed. Lastly, the model is then applied to a real scenario created with a Hele–Shaw cell measurement. The benchmark as well as the comparison with the measurement results show the high potential of this technique.


Fundamental problem
Fiber-reinforced polymers (Voß and Haupt, 2021), highstrength steel , and modern multi-material lightweight construction (Banea et al., 2018) are leading to a steady increase in the use of adhesive bonding technology. On a smaller scale, smart products, equipped with sensors and control units in exposed positions, lead to a significant increase in requirements for sealing electronic components against environmental influences. Both bonding and sealing usually begin with the application of an initially viscous fluid adhesive or sealant, which hardens into a solid in the subsequent curing. The final position of the adhesive or sealant is decisive for the quality of the end product, as transfer of forces requires adequately filled gaps, as designed. Air pockets, or underfilling of the gap, will result in a reduction of load transfer capability and can severely weaken the long-term resistance of the bond due to the ingress of media. Bubbles or interruptions in a sealing bead can lead to the migration of moisture into the component to be protected, where it can cause failure of the overall product. The targeted simulation-based optimization of the adhesive bead position with the aim of completely filling the adhesive gap while using a minimum amount of material is the core concept of an ongoing project of the authors.
The correct position of an adhesive or sealant is not only determined by the accuracy of the adhesive application process, but also by the subsequent step during which the second component is placed on top, creating a squeeze flow between the parts to be joined in the narrowing gap. The current computation of squeeze flows is mostly based on the 3D Navier-Stokes equations. These models solve the flow problem in all its complexity, although only a small subset of the result variables is relevant to the problem. The geometric Vol. 5, No. 1, 2023, 29-52 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-021-0125-3 conditions of an adhesive joint, often a few tenths of a millimeter in thickness while sometimes several meters in length, lead to hefty aspect ratios that results in numeric instabilities and expansive computation time.
In previous studies (Müller et al., 2019), the authors presented a much more efficient model for squeeze flow simulations. The model addresses the aspect ratio and considers only the relevant variables of a squeeze flow in the narrow gap. Originally designed for tribological systems (Müller et al., 2013), it has been successfully reformulated, implemented, and benchmarked for CFD simulations of elementary compression processes in bonding . The model is the basis of the extension presented in this paper to include the influence of trapped air bubbles and pockets.

Significance of air entrapments in adhesives
Defects in the adhesive gap caused by air pockets can have a significant negative effect on the strength of the bonded joint. Various publications deal with the formation mechanisms of air pockets (Guyott et al., 1986;Chester and Roberts, 1989;Davis and Bond, 1999), and conclude that air bubbles may already be present in the adhesive container or may be formed during application or squeezing, trapping areas with air. The effects of air pockets on the strength of adhesively bonded wind turbine blades have been studied in detail in Petryna et al. (2014). Based on experiments and mechanical simulations, it is shown that air pockets generate stress peaks under load, which lead to earlier failure compared to the ideally filled adhesive gap (Petryna et al., 2014;Drieschner et al., 2018). Similar findings by other research groups are presented in the studies published in Zarouchas et al. (2012), Yang et al. (2013), and Ji and Han (2014).
Several publications deal with the influence of defects in bonded joints on their mechanical properties. Here, the focus was mainly on joints made of timber (Grunwald et al., 2014), steel (Karachalios et al., 2013), GFRP (Ascione, 2016), or CFRP (Ribeiro et al., 2016), as well as on improper adhesive mixing ratios (Adams and Cawley, 1988). These works are limited to single and double lap joints. Sengab and Talreja (2016) numerically investigated the influence of air entrapments in bonded joints and found that it is not only the size, but also the shape of the air bubbles that has the greatest influence on crack initiation and growth, with round air pockets less critical here than flat voids.
In the manufacturing of electronic boards, thermally conductive adhesives (also known as underfillers) are injected into a narrow gap (a few micrometers) below a chip to fix the electronic component, dissipate heat, and protect it from corrosion. In such applications, it is essential that no air pockets remain in the gap, as these can severely limit functionality and durability in the long term (Khor et al., 2010). The current electric vehicles use so-called gap fillers to dissipate heat from battery modules into the vehicle underbody (Frauenhofer et al., 2020a(Frauenhofer et al., , 2020b. Here, too, the pressing of these materials in production and the homogeneity of the distribution of the gap fillers play an important role. Experimental and numerical studies have been carried out on this problem (Han et al., 1996;Khor et al., 2010;Chang and Young, 2014), which showed that implementation with model strategies based on similar assumptions used in this work are very effective. It was in particular shown that if the gap height is significantly smaller than the gap width, and that the solid elements within the gap are cylindrical rather than spherical, neglecting pressure gradients and fluxes in the height direction in the simulation still allows to model the experiments very well (Nguyen et al., 1999). All the previous works were performed without accounting for air pockets that may remain in the system even during adhesive injections. This flow behavior, referred to as "Hele-Shaw flow" in many sources, has also been studied and modelled for non-Newtonian fluids (Wan et al., 2007).

Types of air bubbles and existing approaches
Two main types of air bubbles exist in adhesive systems: dispersed air bubbles ( Fig. 1(a), in the following called "air bubbles") and, secondly, areas in which the air forms a large continuous area ( Fig. 1(b)) which are called "air pockets" in this work.
Dispersed air bubbles are present in the form of small bubbles and no large-area phase boundaries. These can occur, for example, on the substrates when a highly viscous adhesive is applied to a very rough surface (Lee, 2013;Ebnesajjad and Landrock, 2014). On the other hand, the air bubbles may already be in the adhesive before application or may enter the adhesive during application (Waugh, 2016). One possibility to describe dispersive two-phase flows is based on the Euler-Euler method. It is then assumed that both phases, the continuous and the disperse phases, form a continuum. For both phases, conservation equations are implemented in the form of mass, momentum, and energy equations. In addition, terms are added which describe the exchange of mass and momentum as soon as the two phases permeate each other. In principle, this method suits modelling separated multiphase flows, but adapted exchange terms cannot be established for all flow types. On the other hand, the suitability for parallel computing is advantageous, whereby a shortening of the computation time is obtained (Spyrou, 2011).
Another modelling approach is the Euler-Lagrange method, in which the continuous phase is described according to the Euler method. The disperse phase, on the other hand, is modelled with the Lagrangian approach. Using Newton's equation of motion, the velocity and location of individual particles are calculated as a function of time (Ferziger et al., 2008). The interaction between the two phases can take place according to the one-way-coupling (influence of the particles on the continuum is negligible), two-way-coupling (particles and continuum interact), and four-way-coupling (additionally particles interact with each other) (Spyrou, 2011). For simulations of adhesive flow using Euler-Euler, Euler-Lagrange, and mixed methods, see Fricke and Vallée (2016) and Landgrebe et al. (2015).
This paper is primarily focused on air pockets as defined above, in which the initial distribution of the adhesive on the part to be joined and the subsequent joining process can lead to their formation. This process can be seen in Fig. 1(b). At the top, the initial distribution of the adhesive before joining can be seen. Then the upper substrate is moved downwards. As soon as it touches the adhesive, an air pocket is created in the center of the bonding joint.
Separated two-phase flows, as they are present for these air pocket types, can be described by splitting them into two single-phase regions. In these regions, the application of conservation equations for single-phase flows is permissible. In these models, the time-varying position of the phase boundary is a fundamental aspect. Due to the fact that a singularity arises in the area of the phase boundary and that abrupt changes of the material parameters are to be expected, a very fine discretization is necessary. Furthermore, the surface tension present at the phase boundary must be taken into account. The methods for describing the phase boundary can be divided into two categories (Ferziger et al., 2008;Spyrou, 2011):  Interface tracking methods: Methods using massless particles as markers of the phase boundary. At the beginning of the simulation, the markers are placed at the known position of the phase boundary, and then their position is tracked. An advantage is the high accuracy and the fact that the number of particles is independent of the lattice discretization. On the other hand, very strong changes of the topology are a disadvantage as they strongly increase the computational effort (Spyrou, 2011).  Interface detection methods: Methods of this type also use markers, but the fluids involved are marked rather than the phase boundary itself. Different approaches for their implementation are reported. In the MAC method, massless particles are added to the fluid, acting as markers as in the interface tracking methods. Because it is necessary to determine the motion of these particles in addition to other equations, a relatively high computational effort is to be expected. The volume-of-fluid (VOF) method, which is used in the CFD simulation presented later in this paper, offers another possibility. Here, a control volume is introduced, which provides information about the position of the phase boundary due to the different phase fractions. Its main advantage is that even complex surfaces can be described efficiently; on the other hand, the indistinct interface between the fluids is a disadvantage. A third approach is the level set method, which describes the surface by a distance function that calculates the smallest distance to the interface at each location in the solution domain. This results in a smooth, linear function in the entire solution domain. This is countered by the disadvantage that no strict conservation of mass is given here, for example, in the VOF method (Ferziger et al., 2008;Li, 2008).

Hele-Shaw flow
A special role of the two-phase flows in the context of this work is played by the Hele-Shaw flow that characterizes the flow in the eponymous cell (Paterson, 1981). In such a cell, as shown in Fig. 2, fluid is located between two flat plates separated by a small gap h.
One type of Hele-Shaw flow with trapped air occurs when there is an opening through which the second fluid of lower viscosity (e.g., air) is applied. This fundamental scenario has been studied by different authors (Baig et al., 2000;Fast and Shelley, 2004). Due to its similarity to this work with regard to the model structure, this flow form will be considered in more detail below. As soon as air is pressed M. Müller,S. Willenbrock,L. Stahl,et al. 32 into the Hele-Shaw cell, an air pocket is formed which displaces the surrounding first fluid. The resulting structure is referred to in literature as Saffmann-Taylor instability, or colloquially as "viscous fingering". Such behavior is always possible when a less viscous medium displaces a medium with higher viscosity.
Modelling approaches for these processes were also discussed in Fast and Shelley (2004). On the one hand, analogies to slow flows in porous media are formulated in the form of Darcy's law (Darcy, 1856;Barton and La Pointe, 2012). If the fluid's velocity is slow enough, then interia can be neglected, and the use of the Laplace-Young equation may be appropriate, in which the interfacial tensions and the radii of curvature of the pockets are of particular importance (Fast and Shelley, 2004).
Further fundamental investigations of the flow of fluids in a Hele-Shaw cell were carried out, for example, in Kondic et al. (1996), in which the non-Newtonian velocity profile over the gap height was considered. Studies on "viscous fingering" during the separation of the plates were published, for example, in Fast and Shelley (2004) and Amar and Bonn (2005). The formation of air pockets in an experiment with very narrowly separated plates in a Hele-Shaw cell was presented using silicone oil and air (Paterson et al., 1995); in this work, evaluation methods related to wetting were also developed. A novel complex model for describing properties of thin films via energetic approaches to wetting can be found in Bertsch et al. (2005), where slip is allowed at the fluid-substrate interface. In addition, an experimental study (Grunwald et al., 2014) investigated the influence of air defects on the load-bearing capacity of bonded wood joints, artificially imitating air pockets in the form of applied Teflon patches. However, fundamental work on gap flows with a focus on simplified assumptions as well as for the description of such joining processes does not exist in the literature.

Generalized partially filled gaps (GPFG) model without air entrapments
The model used, and extended, in the following is based on approaches from tribology; Reynolds derived an equation as early as 1886 (Kondic et al., 1996) to describe gap flows very effectively. In particular, he made two assumptions that made the derivation much simpler: (1) There are no significant pressure gradients perpendicular to the gap plane, which cancels the need for discretization in this direction.
(2) There are only two major momentum components in the Navier-Stokes equation, namely those from the pressure gradients in the gap plane direction, and the viscous resistance due to shear perpendicular to the gap plane. The validity of both assumptions is justified with narrower gaps, if related to the width of the topography, and the viscosity of the fluid. These approaches were taken up by the authors of this paper, and used for fundamental studies on tribological systems with partially filled gaps model (abbreviated as PFGM) (Müller et al., 2013). The approach of separate treatment of the Navier-Stokes equation, and the continuity equation applied therein offer great flexibility for subsequent model extensions (Müller et al., 2014), especially for the problem considered in the present paper. For the application to contact mechanic problems, deformations in the contact bodies were taken into account (Müller et al., 2017). Primarily, in the field of tribology, the PFGM studies have been demonstrated the influence of the amount of air in the gap (Müller et al., 2018a) on the resulting pressures and flows (Müller et al., 2020). Aforesaid simulation techniques have recently been applied to the description of squeeze flows during compression in adhesive systems (Müller et al., 2018b;Kaufmann et al., 2021), where it has been shown that the essential assumptions from the Reynolds equation (and thus also from the authors' model) are fulfilled for the investigation of squeeze flows of adhesives. In a sequel study, non-Newtonian fluids were included in the model . So far, however, trapped air pockets have not been taken into account in this model, which will be presented below. For this purpose, modifications have to be made at different parts of the model.
Before diving deeper into details, directions and coordinates are shown in Fig. 3. In a nutshell, x and y represent the fixed Cartesian coordinates in the gap plane, and z indicates the direction perpendicular to the gap plane. For the following, it is also always assumed that the z-coordinate at the bottom of the gap, thus the lower substrate, is zero. The height of the gap is denoted by h, so any location on the top surface has the z-coordinate h assuming smooth surfaces.
The gap plane directions x and y are discretized, and the assignment to the location is realized by the indices i and j. At each cell (i, j) or location (x, y), different gap heights h and fluid quantities V can be considered, as illustrated in Fig. 4. Figure 5 shows the state variables of the model and their mutual interactions. The variables to be calculated are the  The influence of air in the original PFGM form can be summarized as follows. It does not generate any resistance when it flows, since its dynamic viscosity is significantly smaller than that of the fluid. Air is not explicitly accounted for in the mass balance, so it is not tracked in its motion (refer to Eq. (8) in Appendix A1). Air does not contribute to the pressure buildup; as, if a cell is not completely filled, it is assumed that ambient pressure is present there, cf. Eq. (9). Latter pressure is usually set to 0, so that the pressure p always represents the gauge pressure.
Adhesives represent highly viscous fluids, which are often shear thinning (Pocius and Dillard, 2002), a fluid property for which the PFGM has already been extended for generalized non-Newtonian fluids (abbreviated as GPFGM ). In the extended GPFGM, a subordinate height discretization is implemented at each cell, and thus the mean velocity is derived from the local pressure gradient via the resulting shear rate and flow velocity profile. For further details, please refer to Müller et al. (2021).

Model extension to incorporate the effect of trapped air
All of the aforesaid approaches correspond to the way of modelling tribological systems with partially filled gaps. In terms of trapped air, however, this approach is obviously not appropriate for the following reason. The assumption that ambient pressure is present in cells that are not fully filled is only physically correct if the flow of air can take place via a physical connection to the environment. If this connection is missing or lost, entrapped air may be compressed, and pressure builds up in corresponding cells. Amendment of the model to account for that is derived in this section. The existing model is extended by accounting for the pressure buildup of domains that are not completely filled with fluid and have no connection to the environment. Consequently, a case distinction is to be implemented whether the air in a partially filled cell is trapped or belongs to the ambient pressure area. Before presenting the details of the model in Section 2.2, it is first shown how the detection of trapped air is implemented in an automated way, which is done in Section 2.1.

Identification of trapped air
To explain the algorithm for the identification of trapped air, that is, pockets of air surrounded by adhesive with no contact to the exterior, the main steps are illustrated by an example below. Figure 6(a) shows the two-dimensional top view of an adhesive pattern. The cells shown in red represent fully filled fluid cells, while those in white contain air. The algorithm starts with considering all boundary cells. It is then determined whether the respective boundary cell is partially or fully filled. If it is partially filled, the cell assigns the information that there is a connection to the system boundary, and the ambient pressure p 0 is present in the cell (marked with an x). Therefore, in the visualization of the algorithm ( Fig. 6(b)), all boundary cells are marked in this manner. In the next step, it is now determined whether all cells neighboring the marked ones are also partially filled. If this is the case, the air of these cells also has a connection to the system boundary, and thus the ambient pressure p 0 . The algorithm subsequently analyzes further their neighbors, until all have been considered, as illustrated in Fig. 6(c). If a cell is neither marked by an x, nor fully filled, it must contain trapped air. This is the algorithm's final state leading to the configuration shown in Fig. 6(d), where all cells with trapped air are marked in blue. This algorithm finds all cells with trapped air depending on the distribution of the filling degrees Θ. Due to the change of the fluid volume during the simulation process, this identification is carried out at any time-step. This algorithm is very efficient regarding computing time. Because the filling level distribution also changes during this numerical process within each step, it must be included in the iteration loop of the Newton-Raphson method.

Pressure in areas of trapped air
The model extension is based on a two-phase flow model. Therefore, both the flow of the fluid and that of the air are considered, and the extended model contains further state variables for the air flow velocities and the volume of air. It is assumed that the pressure in a mixed cell is homogeneous and equal for both phases. The summarized set of state variables is shown in Fig. 7(left).
In order to distinguish the variables associated with air and fluid, the former is labeled with the index G (as gas) while the latter is given the index F (as fluid). It should be emphasized at this point that the volume state variables always represent uncompressed volumes, so that the correct masses are actually accounted for in the continuity equation. The volume balance in the compressed state is of particular importance here, and thus Fig. 7(right) defines the volume terms used as follows. The compressed volumes are additionally assigned the index C; the volumes in the uncompressed state are thus composed of the compressed volumes and the respective difference volumes, as shown in Eqs. (1) and (2).
Consequently, the total gap volume consists, as given by Eq.
(3), of the two compressed volume fractions.
Since the model already works with compressibilities, the influence of compressed air can directly be taken into account. The air is implemented as an ideal gas (Clapeyron, 1834). From the viewpoint of thermodynamics, compression and heat generation is a polytropic process; if the compression is very slow, the process is isothermal, as generated heat has enough time to dissipate. Numerically this is the simpler case for modelling air, as detailed in Appendix A2. If compression is performed so fast that local temperature rises, then the process is adiabatic. A generalized approach for a polytropic process is the subject of Appendix A3. Practical cases will be somewhere between both idealisations. The general form of the new model can be found in Appendix A4.
The final equation for the pressure in a cell (i, j) is summarized in Eq. (4), where n is the polytropic index. It is 1 for an isothermal process, approximately 1.4 for air within an isentropic process and can take different values in between for a polytropic process. , i j Θ is the local filling level, which is the ratio of fluid height to gap height. E is a penalty factor that represents the compressibility of the fluid, and p 0 is the ambient pressure (more details are documented in Appendix A3).

Momentum and mass balances in areas of trapped air
In addition to the fluid flows generated from the pressures from Eq. (7) in Appendix A1, the flow of the air must also be taken into account for systems with trapped air. Air cells with a connection to the environment are skipped as their corresponding flow is not of interest. The flow velocities, which describe the air flow, can now also be described via the simplified Navier-Stokes equations accounting for the viscous term and the pressure term (momentum balance). The terms are then included in the continuity equation (mass balance) (Müller et al., 2014). With the viscosity of the air, the following equations result for the description of the air flow velocity (see Eq. (5)). The viscosity of the air was set to Looking at the equations, it can be seen that an air flow velocity is calculated whenever there is a pressure gradient between two cells. The same applies to the velocities of the fluid flow. However, the ratio between the fluid and air flow is unknown and can not be determined a priori. Therefore, either further equations are required for the phase boundaries or case distinctions must be made to define which flow occurs. The aim of the model presented here is not to represent all processes with the greatest precision, but to use a pragmatic approach to track the influence of trapped air on the system behavior and, in particular, on the adhesive flow. With this in mind, case distinctions are used here.
This case distinction is based on the states of the discrete cells. Figure 8 depicts an exemplary, one-dimensional scheme, including trapped air, free air, and fully or partially filled fluid cells.
In this and all following illustrations, the fluid is shown in red, trapped air (if applicable) in blue, and air with a connection to the environment in white ("free air phase"). For Fig. 8, the air (blue color) in the middle (cells 3-6) is trapped, because the cells do not have a connection to the environment. However, the air in cells 1 and 8-10 is not trapped and will not create a pressure buildup.
For the flow between each of the cell, a variety of different rules apply. Figure 9 shows a corresponding overview.
For each case, two adjacent cells are shown (onedimensional case), and the flow across the boundary of the cells is considered. Some examples referring to the scenario in Fig. 8 are: case 2 for the boundary between cells 9 & 10; case 4 for cells 2 & 3 or cells 6 & 7; case 5 for cells 4 & 5.
Below each schematically illustrated case, the labels with green check marks, and red wrong cross symbols, represent the medium flowing in the respective case, and the arrows represent the associated directions.
Case 1 represents two cells fully filled with fluid. Therefore, only fluid flow velocities are to be calculated between these cells. In case 2, a partially filled cell and a fully filled cell are shown. Since the partially filled cell does not contain any trapped air, only fluid flow is considered. Case 3, two adjacent cells with trapped air, leads to air flow.
While cases 1-3 do not need any case distinctions, cases 4-8 show scenarios where both phases could potentially flow. To limit the model's complexity, a pragmatic and efficient approach was chosen here taking into account the greatly differing viscosities and compressibilities of air and adhesive. During accompanying experiments it was found that the phase boundary moved over the entire gap height (vertical phase boundary). Flow patterns where the air pocket flowed above the fluid were not observed. This led to the general assumption that the phase boundary only spreads across one discrete cell (between the entrapped air domain and the fluid domain there is maximum one mixed cell). From this general rule, more detailed case distinctions were derived for the specific cases. These distinctions yielded good results for the investigated systems of this work. However, different fluids may show different behaviour.
In case 4, the flow between a mixed cell and a fluid-filled cell is considered. With the mentioned simplification in mind, in this case, only the fluid flows. On the one hand, this mechanism prevents the air from connecting to the pressureless edge region through a very thin gap, and on the other hand, it prevents an accumulation of fluid volume "under" an air pocket. For the benchmark in Section 3 and also the practical scenarios presented in Section 4, this represents well the measured behaviour.
Case 5 concerns two cells, both containing trapped air. In addition, the left cell has a fraction of fluid volume. For this situation, it is assumed in the model that only the air can flow between the cells. This can be justified by the fact that the air flows faster due to its significantly lower viscosity and is, therefore, more likely to provide pressure equalization than the fluid. For the same reason, also in case 6, in which two partially filled cells with trapped air can be seen, only air flow and no fluid flow occurs.
For cases 7 and 8, a fully filled fluid cell neighbours a cell with trapped air. In case 7, there is a pressure gradient from the right to the left cell. In such a case, a flow of air over the boundary of the cell is calculated. However, if the pressure gradient should be in the other direction, as is the case 8, a fluid flow velocity must be calculated.
With the required velocities now defined, the volume balance for the air can be formulated analogously to the fluid volume balance equation in Eq. (8) (in Appendix A1), to yield Eq. (6).
Overall, the chosen model is a pragmatic approximation of the squeezing flow with entrapped air. The assumptions used are based on experimental observations such as during adhesive compression processes the air pockets do not flow infinitely fast, nor do they regularly flow in a mixture with the fluid. The air flow is thus decoupled from the fluid flow. A trapped air pocket can thus continuously change shape, since a flow of air into previously fully fluid cells is possible.

Benchmark test
At first, the model will be tested with a benchmark simulation in order to verify the individual assumptions introduced in Section 2. For this purpose, the PFGM solution is compared to a CFD calculation for the case of a squeeze flow under the influence of trapped air.

Investigated scenario
The scenario considered, shown in Fig. 10, is a two-dimensional problem in terms of CFD, and a one-dimensional problem for PFGM. It models two rectangular adhesive beads with an initial width of 5 mm, which are infinitely extended in depth, and between which there is an area of trapped air At the left and right system boundaries, there is ambient pressure. The gap height 0 h at the beginning is 2 mm. At the beginning of the simulation, there is ambient pressure in both the trapped and free air. The lower wall is fixed, and the upper wall moves downwards at the constant speed of 9 mm/s to a final gap height of 0.7567 mm. This final value is obtained from the last convergent step of the CFD calculation with automated time-stepping before the fluid leaves the domain. The total process time T is therefore Δ / 1.2433 mm (9 mm/s) 0.13817 s T h w = =¸= . The criteria for the subsequent comparison are the pressure distribution in the horizontal (or x-) direction, and the mean pressure over time. In addition, the volume of the trapped air pocket, as its average width, is examined.

CFD model simulations
The benchmark scenario is firstly computed using the commercial CFD software Ansys 2020 R1 Fluent. Both walls are set to be adiabatic with the heat flux of zero. A uniform structured mesh is applied for the whole domain. The mesh size (40 elements vertically, 600 horizontally) was chosen after a sensitivity analysis to ensure convergence and accuracy of the results. The volume of fluid (VOF) modelling strategy is implemented for this multiphase problem. The energy equation is also activated so that the thermodynamics of the system is taken into account. The extensive details about the dynamic layering mesh, and all settings for the solver, were fully described in Müller et al. (2018b). The time-step size is automatically adjusted through the adaptive time- Towards the efficient modelling of trapped air pockets during squeeze flow 37 stepping algorithm. The simulation requires approximately 10,000 time steps, which represented 7 hours on a standard workstation. The results are summarised in Fig. 11, which shows the motion of the fluid in the gap, and the resulting pressure, at selected time steps. The pressure indicates the gauge pressure. Figure 11 shows that at the very beginning (t = 0.02188 s), the fluid moves outwards, but also slightly in the direction of the air pocket, causing it to decrease not only in height but also in width. This is accompanied by a pressure buildup in the air pocket. At the second instant (t = 0.08888 s), the air is compressed to such an extent that the fluid flows exclusively outward. This trend continues, but as the gap height decreases, the resistance of the fluid against shearing also increases, which is why the air pocket becomes thinner again (t = 0.10005 s), which continues (t = 0.13817 s). The fluid distributions also show that the fluid takes on a more rounded shape at the front of the outer regions. The fact that there is no fluid in the outer regions directly at the wall is mainly due to the no-slip boundary condition. This effect, already described in detail in Müller et al. (2018b), is not likely to represent the true processes (in which the wall would be wetted), but rather represents an inaccuracy of this CFD model; consequences resulting thereof will be discussed in Section 3.4.
No significant pressure gradient occurs over the height (vertical direction) in the fluid, which backs the corresponding assumptions of the PFGM. Furthermore, as expected, there is no pressure in the free air, while built up in an approximately constant distribution inside the air pocket. This will also be discussed in more detail in the upcoming Sections 3.3, 3.4, and 4.2.2.
In the CFD model, the thermodynamics of the system is covered by the energy equation. To enable a comparison with the PGFM, the polytropic index is analyzed in more detail. To do so, the average pressure in the compressed air, and the compressed air volume, are extracted from the model at each point in time. From this, a pressure-volume diagram can be generated as shown in Fig. 12 (see red curve). With the initial volume and the ambient pressure existing at that time, corresponding comparison curves for the isothermal ( 1 n = ) and isentropic case ( 1.4 n = ) can be obtained; see also Eq. (24) in Appendix A3. The corresponding polytropic index was fitted and set at 1.23 n = , which is consequently used for the comparability in the PFGM. The air temperature rose from 300 to 384 K.

PFGM calculations
In the PFGM, the domain was discretized in 1D with 300 cells (in the x-direction). The time-stepping is adaptive, and 437 time steps were required for the total time; on a comparable workstation it only took 15 min, thus around 1/28 of the CFD computation time. The polytropic index 1.23 n = was set. As in all previous PFGM applications, pressure and velocity changes between successive iteration steps, and mass conservation, are formulated as convergence criteria; very robust convergence was observed in all cases. Depending on the fluid distribution, only 10-20 steps of the Newton-Raphson method were needed for each time step. Figure 13 shows the pressure curves and the fluid distributions in the first phase of the process.

Time (s)
Fluid/air distribution Pressure distribution 0.02188  The pressure curve over the longitudinal coordinate (x) at varying instants of time is shown in Fig. 13(a), while Fig. 13(b) depicts the width of the air gap versus time. This is calculated for the region in the center with trapped air as follows. The algorithm detailed in Section 2.3 results in a scenario in which a cell filled with compressed air and fluid must have a neighbored cell completely filled with air and a neighbored cell completely filled with fluid. The air gap width is thus calculated by the sum of connected cells completely filled with compressed air and the respective air fraction of the partially filled cells (that represent the phase boundary). This interpretation represents a phase boundary that is orientated in the vertical (z-) direction. Figure 13(a) shows that at the beginning (t = 0.00033 s, blue line), the pressure inside is already higher than zero due to the compression, but it is still significantly smaller than the fluid pressure (as can be seen from the peaks on the left and right). Consequently, the fluid flows in the direction of the air pocket according to case 8 and later to case 4 as defined by Fig. 8, which is also indicated by the reduction of the air pocket width, see Fig. 13(b). With increasing time, the pressure in the entire domain increases. Due to the further compression of the air pocket, the pressure in the air increases faster than in the fluid. From about t = 0.00333 s, purple line, the location of the maximum pressure is already in the air pocket; as expected, it is constant. By increasing the pressure in the air pocket, it occurs that, in accordance to case 5 and later also case 7, the air flows to the outside, i.e., the air pocket widens. This behavior can also be very well noticed in Fig. 13(b), where from approximately t = 0.0032 s the curve rises again. The algorithm, with its case distinctions, thus works plausibly for both the reduction and the increase of the air pocket width.
With further compression, the pressure in the air pocket increases considerably, as shown in Fig. 14. In the process, the fluid now also flows more intensively to the outside, so that further fluid pressure domains are created there, see Fig. 14(a). As the gap becomes narrower, the resistance of the fluid against shear increases, so that the pressure in the fluid is again greater than in the compressed air (starting between t = 0.09967 s and t = 0.1155 s), which leads to a renewed flow of the fluid in the direction of the air pocket and causes it to become thinner again. This behavior can be seen accordingly in Fig. 14(b) from a time of approximately t = 0.9 s. In summary, it can be concluded that the PFGM works as intended in this benchmark. The extent to which these dynamics are consistent with the CFD calculations will be discussed in Section 3.4.

Comparison between CFD results and PFGM results
Pressure and air pocket width are again used to compare the results obtained by CFD to PFGM detailed above. Figure 15 shows the average pressure over the entire domain (including the zones with free air, Fig. 15(a)) and the width of the air pocket ( Fig. 15(b)), both plotted against the gap height. The values can also be interpreted directly over time considering the initial height of 2 mm and the constant vertical velocity of 0.9 mm/s. The width of the trapped air pocket was calculated in the PFGM as described in Section 3.3, in the CFD using the air volume present at the instant (as a summation over the air fractions in the center area) divided by the specified gap height. The comparison of the pressure ( Fig. 15(a)) shows identical, highly progressively increasing trends, where PFGM values lie above those obtained by CFD in the entire interval. The deviation, expressed as relative difference, increases as the gap height decreases, and by the final step it amounts to 20%. The comparison of the pocket width also shows very similar results. In particular, it is clear that all trends are identical, i.e., the slight decrease in pocket width at the beginning, then an increase and finally the decrease again. Again, as the gap height decreases, the difference increases, but the relative difference of about 10% towards the end is still relatively small. The larger air pocket calculated in the CFD correlates clearly with the smaller pressure. This leads to the conclusion that the case distinction (Section 2.3) for the phase flow is valid.
To identify the cause of these differences, the curves of the fluid distribution and the pressures over the x-axis at different instants are examined (see Fig. 16).
On the one hand, this concerns the distribution of the fluid in the gap represented by the value for θ (Fig. 16(a)). In the CFD case, an average value of the fluid fraction over the height (z-direction) was calculated for the respective x-position. For the pressures (Fig. 16(b)), also the mean value over the z-axis was determined and plotted for the CFD calculation. Since the time steps differ between CFD and PFGM, the CFD instants were used as a basis, and the curves from the PFGM results were interpolated between the corresponding adjacent instants.
In the curves, for the filling ratio θ , it can be seen that the central region of the trapped air pocket agrees well in both calculations. Larger discrepancies are found in the outer regions, especially in the regions where no fluid was initially present. As explained in Section 3.2, in particular Fig. 11, this is a consequence of the no-slip condition on the walls resulting in no fluid ever reaching the walls. Since, as mentioned earlier, this behavior does not correspond to physical experience, the solutions from the PFGM should be much more accurate in these regions (as these regions are also fully filled). This effect has also already been studied and interpreted in detail in Müller et al. (2018b). The pressure curves depict a very good agreement, with slightly increasing ).

Fig. 15
Comparison between CFD and PFGM for the compression process.
deviations with further squeezing. Due to the lack of wall contact of the fluid in the outer areas, the pressure here does not increase as much from the outside to the center. As a result, the discrepancy between CFD and PFGM also increases with increasing distance from the outer boundaries. Generally, however, it can be stated that the PFGM describes this benchmark very well and that the deviations from a much more computationally expensive CFD calculation are relatively small. It can even be assumed that, depending on the scenario, the PFGM partially generates more accurate results than the CFD, as the CFD calculation used here only weakly describes certain physical effects such as the flow near the wall. The comparison of a PFGM calculation with measurement is discussed in Section 4.

Compression and depressurization of trapped air pockets
In real systems, during pressure buildup, the air will usually move in the direction of a weak point (i.e., an area where  the surrounding fluid width is small). When it reaches the boundary with the ambient air, a pressure release will occur in which the portion of compressed air escapes from the air pocket and only uncompressed air remains in the pocket. The model will be tested below to demonstrate that scenarios in which an air pocket opens and closes again accordingly (Section 4.1), and an experimental test is compared with a simulation (Section 4.2).

Simulation for opening and re-closing of air pockets
As already detailed in Section 2.1, identification of whether trapped air is still trapped is performed after each time step. In this context, the cells that are in contact with the ambient pressure boundary can be specifically detected. If this is the case, the difference between uncompressed and compressed air volume G , Δ i j V is removed from the corresponding cells. This process imitates an infinitely rapid pressure release that occurs as soon as there is no cell with fully filled fluid between the air pocket and the ambient air.
The decompression is shown for an example simulation in the sequence of Fig. 17. In these figures, the cell state is shown in a two-dimensional view in the top row and the corresponding filling levels (in terms of the uncompressed volume) are shown in the bottom row. Consequently, the air that exceeds the fluid level corresponds to the difference between the uncompressed and compressed air volume.
The initial configuration at time step 0 (not shown here) is an adhesive pattern in the form of a ring made thinner at one spot. In Fig. 17(a), the cell state at time step 230 is shown. This instant corresponds to the step right before the depressurization occurs. It can be seen that between the air pocket (blue) and the ambient pressure cells (white), there is at least one cell completely filled with fluid (red). During the next step (time step 231, Fig. 17(b)), there is now a location where the air pocket comes into contact with the ambient pressure area. As a result, the corresponding compressed air volume is released at this time step. The pressure release and its corresponding flow is not explicitly computed. It happens much faster than the model's time step. The formerly trapped air phase (compare to Fig. 8) changes its state to free air. Consequently, in time step 232 (Fig. 17(c)), no more compressed (and trapped) air can be detected in the domain.
As a result of the depressurization, there is no longer any overpressure inside the pocket, and the fluid subsequently flows in the direction of the original pocket. In this case, this leads to the formation of a new area with trapped air within the following time steps. However, the air is obviously not as compressed as before the depressurization. This process is illustrated in Fig. 18. At time step 244 ( Fig. 18(a)), there is still a connection between the inner and outer areas. This channel is closed by the fluid movement during time step 245, see Fig. 18(b). Consequently, the algorithm from Section 2.1 detects again trapped air (blue), thus (according to the case distinction in Fig. 9) changing the situation at the boundary from case 2 to case 4. In the following time step 246 (Fig. 18(c)), the vertical motion of the upper surface compresses the trapped air again (seen by the small volume exceed of air volume compared to the fluid volume). Consequently, these real processes can be represented well with the aid of the model. Associated quantities for pressures during these processes are presented in Section 4.2.2.

Measurement procedure and results
In order to validate the motion and the influence of trapped air pockets by measurement, an experimental setup is used that is based upon a Hele-Shaw cell represented in Fig. 19(a), which is mounted in a universal testing machine ( Fig. 19(b)).
The Hele-Shaw cell consists of two glass panes that can be translated in parallel to each other, with the distance controlled by the universal testing machine (UTS). This setup is described in more detail in Kaufmann et al. (2021). For the tests presented herein, a high-viscous silicone oil (Elbesil B 2.000.000, QUAX GmbH, Germany) was chosen, as the fluid well represents a conventional adhesive without time-dependent effects related to curing. This high-viscous fluid has been rheologically investigated in advance with a rotation rheometer and a capillary rheometer, so that the shear dependence of the viscosity was determined in addition to the viscosity . The rheological characterization of the fluid is illustrated in Fig. 20.
For shear rates smaller than 1 s −1 , the silicone oil exhibits nearly Newtonian flow behaviour. The corresponding viscosity amounts 1744 Pa•s. In this paper, the model is initially applied only to Newtonian fluids. For the system investigated in the following, it was subsequently found that the maximum shear rates of 10 s −1 can occur at the boundaries in the final position; on average over the entire domain, the effective shear rate is about 3 s −1 . For this shear rate, the viscosity according to Fig. 20 is about 1500 Pa•s, i.e., about 15% less than the viscosity for very small shear rates. This results in about a 15% overestimation of the load in the model compared to the experiment for the final steps of the simulation. This is considered acceptable in light of the  high complexity of the system. Further discussion on this matter is given at the end of this section.
Before the measurement, the Hele-Shaw cell's glass panes were cleaned with isopropanol. The fluid was then applied to the lower glass pane using commercially available disposable syringes. The initial application pattern was a ring-like pattern so that air entrapments were generated. The amount of fluid applied was determined by weighing. The upper glass pane is moved downwards in conjunction with the support plate and the movable crosshead. A camera mounted above the glass panes records the flow behavior in the gap plane directions. A force sensor is also attached to the support plate to measure the normal force generated by the viscous resistance of the fluid. In addition to the force, the distance between the pane surfaces, so the gap height, was measured. The joining process was then completed starting from an initial gap height of h = 5 mm. The vertical velocity was 0.1 mm/s. The compression process was finished at a final gap height of 0.4 mm.
The fact that the initial application was carried out manually means that the test can only be reproduced to a limited extent. In comparison with the simulation, the focus is on how a pocket moves from an initial configuration with air pockets, and whether the sudden jump in force during pressure compensation can be reproduced. The experiments must therefore only allow a qualitative description of the flow behavior. Exact reproducibility of the application patterns and an exact repeatability of the experiments was therefore not mandatory in these investigations. The further explanations therefore describe frequently occurring phenomena by means of an example. Figure 21 shows a sequence of images documenting the fluid distribution during the compression process. The image at time t0 shows the moment when the upper glass pane first touches the fluid (to get a better impression of the time intervals, the corresponding value here has been set to t 0 = 0). In the center of the fluid air is trapped; additionally, it can also be seen that the fluid contains a low amount of small dispersed air bubbles.
At the next point in time t 1 , fluid has already spread outwards as a result of the squeezing process, and the air has already undergone compression. This is followed by a process where air already spreads (at t = t 2 ). The air inside is thus compressed to such an extent that the resulting pressure causes the fluid to be displaced. The image at time t 3 shows the instant immediately before the pocket depressurization which is accompanied by a sound of an air pocket bursting. It is also notable that the air pocket has a slight contraction, which causes the air pocket to split into two separate air pockets (at t = t 4 ). In addition, it can be seen here that the randomly dispersed air bubbles are now barely visible in the fluid. When the final gap height is reached (at t = t 5 ), most of the air has left the fluid. The remaining air pocket is very close to the boundary of the fluid domain.
The load-gap diagram is shown in Fig. 22. A significant  force is generated from a gap height of about 2 mm on, which then increases progressively with decreasing gap height. This phenomenon has been studied, and presented, in detail in previous publications (Müller et al., 2018bKaufmann et al., 2021). The depressurization of the air pocket is characterized by a sharp and sudden force drop of approximately 2 N; it occurs at a gap height of approximately 5.4 × 10 −4 m. After the drop that resulted from the pressure release, force continues to rise again until the final gap height is reached. The occurrence of such a force drop associated with a pressure release of a pocket was well reproducible in further tests.

Simulation with the PFGM and comparison with the measurement
For the presented experiment to be simulated with the PFGM, the initial distribution must be rendered. For this purpose, the video footage is evaluated at the time of the first contact between the glass pane and the fluid (instant t = t0). A grid of 40 × 40 cells is chosen as discretization and optically overlaid on the fluid distribution, as shown in Fig. 23(a). Then, the distribution is rotated by 45° in order to design the fluid flow in the direction of the weak point parallel to one of the two coordinate directions. The development of pressure through the trapped air and the associated system properties are illustrated in Fig. 24, which exemplarily represents the situation during the compression phase of the air at the gap height of 1 mm (which is between t1 and t 2 ).
In Fig. 24, the pressure distributions in the fluid (top left) and the air (bottom middle) are shown separately. The bottom left is shown the velocity of the fluid. At the top middle, the current position of the upper disk is drawn in weak yellow, as well as the filling level of fluid in red and trapped air (in blue), both in terms of the uncompressed state. In the upper right is shown the associated binary membership, whether a cell is full of fluid (red), partially or completely filled with trapped air (blue), or partially or completely filled with untrapped air (white). The bottom right image shows the local filling ratio of the fluid.
During the compression process documented in Fig. 24, the simulation shows a pressure buildup in the air pocket which exists because it is compressed to almost half its volume (see top middle). As to be expected, the pressure inside the pocket is almost constant (bottom middle). The corresponding value meets the fluid pressure (top left) at the boundary between the phases. This situation also leads in the model to the case of air displacing the fluid (see Fig. 9, case 7).
The higher the compression, the greater the resulting pressure, so that the air pocket consequently spreads further in the direction of the weak point (to the right). This change can be clearly seen when the initial fluid distribution (Fig. 23) is compared with that shown in Fig. 24 (bottom right). The propagation in this direction occurs because the lowest backpressure of the fluid is present here. This process continues until the air pocket reaches the ambient pressure boundary and thus a depressurization occurs. The simulation results at that time step can be seen in Fig. 25. After depressurization, the pressure in the fluid drops, and the empty cells inside the ring fill with fluid. Subsequently, the fluid flow at the opening causes the air pocket to close again and, consequently, the pressure inside the air pocket and in the fluid to increase again (see Fig. 26). Figure 27 shows the comparison for the fluid distribution, at the moments t0, t 1 , t 2 , and t 3 . It can be seen that the air pocket in the simulation at time t 1 is already slightly more spread out in its elongated shape than in the experiment. Such a slight deviation in the shape of the air pocket can also be seen at time t 2 . In the simulation, the round shape seems to change into a rather elongated channel in the course of the compression process, while in the experiment it appears more  compact. As will be shown below, depressurization occurs somewhat earlier in the simulation than in the measurement. At time t 3 (which is after the first depressurization of the simulation and before the first depressurization of the measurement), the size of the air pocket is therefore significantly smaller in the simulation than in the measurement.
This sequence of images shows that the direction of movement of the air pocket is very similar, and that the depressurization occurs in each case at the location of the initial fluid weak point that corresponds to the area with the smallest fluid ring width. The differences in the evolution of the shape of the air pocket and the moment of depressurization can be due to model inaccuracies, or deviations in the initial fluid distributions.
Finally, the calculated normal force curve of the PFGM is compared with the force curve to that of the experiment in Fig. 28. As shown previously, the pressure drop, or deflation, in the PFGM occurs at a slightly larger gap height than in the test. However, this difference of about 0.02 mm is relatively small considering the complexity of the system.  The absolute values of the normal force, and its pressure drop, by about 3 N are similar to the measured of 2 N. Considering that, as mentioned above, the viscosity in the model is about 15% higher than the real system for the last time interval, the existing deviation between the pressure drops would become even smaller if shear thinning is taken into account.

Conclusions and outlook
This paper presents the prediction of flow processes of adhesives during compression in the presence of air pockets. The model is based on the very efficient and established partially filled gaps model (PFGM). The great advantage of this approach is the reduction of the system by one dimension, since the direction perpendicular to the gap plane does not have to be discretized. An extension was developed and implemented to model the compression of air and incorporating the effect of trapped air pockets. The presented approaches also allow the consideration of a polytropic compression, i.e., systems with polytropic index unequal to 1.
The model was first compared with a CFD calculation used as a benchmark. This showed very good agreements with regard to the contraction and expansion of the trapped air pocket. Also, the resulting pressure curves with alternating maximum pressure in the air pocket and in the fluid have been reproduced very well by the PFGM. The slight discrepancies between CFD and PFGM are acceptable and, as shown, the PFGM solutions are probably more accurate.
To simulate the movement of an air pocket in a real system, an algorithm was developed enabling compressed air in air pockets to escape by pressure equalisation when the pocket connects to the surrounding area. After the air pocket has been depressurized, the associated air region has not yet disappeared, and a new pocket may also be created at this location under certain circumstances.
Finally, the extended model was compared with a measurement in which a fluid ring was initially created. This initial distribution was also used in the simulation. The motion of the pocket in the direction of the thinnest fluid point, as it occurs in the experiment, has been demonstrated in the simulation. In addition, the force jumps during the depressurization of the pocket are similar in magnitude.
Overall, it can be stated that the PFGM can be used to model the process of pocket compression, pocket motion, and the resulting force curves. The viscosity of the fluid has a major influence on the processes explained in this paper. As the viscosity increases, the associated pressure usually increases as well, since the fluid poses a greater resistance to the flow and thus increases the force needed to bring the glass panes or substrates closer together. On the other hand, the dynamics related to air pockets is directly influenced by the viscosity for the following reason. At the boundary between the air pocket and the fluid, there exists an equilibrium between the pressure due to the compressed air and the pressure due to the viscous flow resistance, leading to larger compressions (and reduced flow) for higher viscosities. This will affect the formation of the entrapped air pocket and its opening and closing during squeezing. These effects will be investigated in an upcoming paper.
The PFGM will thereafter be used for systematic studies of squeeze flows during the compression of adhesives with complex application patterns, in order to identify an optimum between applied adhesive quantity and a minimum of remaining air pockets. Furthermore, this technique will not only be used for squeeze flows but also for injection-driven bonding.
For both problems (squeeze flow and injection), dispersed air bubbles can in principle also be treated. Since these move predominantly with the fluid, their influence can be considered primarily with a modification of the fluid properties in the model, without the need for further extensions as for air pockets. Corresponding modifications have already been tested and will be the subject of future studies.
Furthermore, the new model can also be used to investigate the influence for different fluid properties like density or non-Newtonian viscosity. Furthermore, with reasonable effort, extensions to the model can be made to include correction parameters and flow factors so that the model can account for surface wettability or rough surfaces. These studies will be published in a future research paper.
The volume in a cell is in turn correlated with the pressure in the cell. In the context of the existing model, the strategy of artificial compressibility is used. For further details, see Nithiarasu (2003), Eq. (9).

( )
If the cell is not completely filled with fluid, the pressure is assumed zero. If a cell is overfilled, i.e., Θ > 1, the pressure is built up, with the amount of overfilling being proportional to the pressure, and the constant E acts as a proportionality factor. This factor is infinite for incompressible fluids and represents the compressibility of the fluid if the latter is to be modelled. In this sense, this pressure-volume relationship acts as a penalty function to a constraint equation.
The resulting equations are a set of implicit nonlinear algebraic equations (Eq. (10)), to be solved with the Newton-Raphson method (Müller et al., 2018b). For an isothermal process with an ideal gas, the product of pressure and volume, pV = nRT, is constant (Clapeyron, 1834). In the case that the uncompressed air volume is taken as a reference, which exists at ambient pressure 0 p , aforesaid ideal gas law can be rewritten to Eq. (11).
First, a cell that is fully filled with air is investigated. Thus, the compressed air volume is equivalent to the gap volume ( G,C , ); Eq. (11) be directly reformulated for the absolute pressure a p as Eq. (12).
However, the present model should basically work with the gauge pressure (which is marked without any additional index), for which thus applies: For a cell completely filled with fluid the differential In the mixed filling case both compressions have to be considered at the same time. The latter two formulas can also be balanced to any other volume, thus also to the respective compressed volume of the fluid/gas, as shown in Eq. (16).
The compression components can be summarized according to Eq. (17).
The proportion of the fluid volume to the total volume can again be described via the filling ratio Θ introduced previously in Eq. (7), leading to Eq. (18).

F , ,
Gap , For simplicity, it is now assumed that due to the significantly higher compressibility of the fluid compared to air, this value is approximately the same as the value for the uncompressed fluid volume, resulting in Eq. (19).
This assumption is only to be questioned for values close to 1 Θ = , which will be discussed in more detail at the end of the derivation. Under all these aforesaid assumptions, the fraction of the compressed air volume can also be given by Eq. (20).
( ) Equation (23) matches the form in Eq. (9), with the difference that here the sum of fluid volume and air volume is used and the proportionality factor is a compression modulus weighted by the fluid and air fractions. At this point, it is worth mentioning again the assumption that for the calculation of the value of Θ , the compressed fluid volume was used instead of the uncompressed volume.
For the critical case of 1 Θ = (fluid only, full filling), the value E would result for the , i j K according to Eq. (22), which means that the compression modulus corresponds to that of the fluid, which is correct. The actual overfilling and the resulting pressure is, after all, calculated with the uncompressed volume, so that Eq. (9) is obtained again for the case of the gap completely filled with fluid.

A3 Pressure equation for the generalized approach
The compression during industrial adhesive bonding is rarely isothermal. Therefore the compression can be modelled with a polytropic process. Here temperature, pressure, and volume are linked via Eq. (24) (van Wylen and Sonntag, 1985).
Here, T stands for the temperature, the index 0 for the initial state, and n for the polytropic index. A polytropic index of 1 is identical to the isothermal case. The value for n is approximately 1.4 n = for air at ambient temperature (Lange and Dean, 1985). Following the same procedure for the derivation of the pressure relation formulated for the PFGM, Eq. (25) is obtained for air-compare to Eqs. (13) and ( As a result, the compression modulus becomes larger and the system therefore stiffer. This factor is 1 / n for very small gauge pressures, i.e., about 0.71 for air, and 0.25 for a pressure of one hundred times the ambient pressure ( , 0 / 100 i j p p = ), so α is within one order of magnitude for realistic scenarios. Since the system of equations to be calculated is implicit anyway, this further dependence of the pressure on the pressure does not lead to any major numerical difficulties. The pressure from the current iteration step is used ( Δ , ) t t j i p + to calculate the compression modulus.