## 1 Introduction  
Thin-walled box girders are affected by a phenomenon known as the shear lag effect, which is the non-uniform stress distribution along flanges caused by shear deformation [1, 2]. This effect increases deflections and peak stresses, causing bridge cracking and collapsing [3]. This topic is well known and widely studied in the literature [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Moreover, in recent decades, many scholars have adopted variational analyses to study the shear lag effect [14, 15]. Although the variational energy method is simple and well applicable, its accuracy depends on the selected shear lag warping displacement functions. Many authors have contributed to this. [PERSON] [14] was the first to adopt the variational energy method to study the shear lag effect and established a quadratic parabolic model to describe the normal stress distributions of rectangular box beams. However, other scholars discovered that the shear lag warping displacement function with a cubic parabola was more suitable for analyzing the shear lag effect of box girders [16, 17, 18]. Combining the above two models, [PERSON] [19] and [PERSON] [20] proposed that the warping displacement function of top plates was a quadratic parabola, while the warping displacement function of cantilever plates was a cubic parabola. In addition, cosine function models were also used to address the shear lag effect problem [21]. An additional term was introduced to the cosine function model to ensure the axial equilibrium of the single-cell box section [22], which improved the accuracy and preciseness of shear lag analysis. [PERSON] [23] then defined two global modification factorsthat prevented the generation of additional axial force and bending moment caused by shear lag warping stresses, successfully capturing the self-balancing characteristic of the system. Additionally, the web's longitudinal displacement was added to satisfy axial equilibrium conditions, and independent shear lag warping displacement functions were employed in different flanges [24; 25]. [PERSON] [26] also considered axial force balance when formulating the shear lag displacement model of the single-cell box girders and pointed out that considering axial force balance yields more precise results. However, although there have been many assumptions regarding displacement function, little attention has been paid to analyzing the shear lag effect in box girders from the perspective of shear flow distribution, particularly in double-cell box girders. [PERSON] [27] focused on single-cell box girders' shear flow distribution law and defined the shear lag warpage function through the shear deformation law. [PERSON] [28; 29] redivided the flange width of multi-cell box girders based on the shear flow distribution before depicting the stress distribution with the cosine function. Nevertheless, the above research reviews could not obtain precisely accurate normal stress distributions, as they failed to reflect stress differences at different webs for the double-cell box section, and the process of calculating the coefficients was complex. Additionally, there is limited research on the shear lag effect in double-cell box girders. The difference in shear lag coefficients for different web plates is not reflected in [30; 31]. The authors of [32] present only finite element simulation results and do not offer specific displacement functions. This paper presents an analytical method based on the shear flow distribution law to study the shear lag effect of thin-walled single- and double-cell box girders. The first step in this method is to determine the box girder's shear flow distribution. Subsequently, a series of novel improved longitudinal displacement functions mathematically expressed as cubic parabolas are established.