Robust Numerical Solution for Solving Elastohydrodynamic Lubrication (EHL) Problems using Total Variation Diminishing (TVD) Approach

TL;DR

This paper introduces a stable TVD scheme with hybrid line splittings for solving complex elastohydrodynamic lubrication problems efficiently, reducing computational complexity and handling high-order discretizations without perturbing robustness.

Contribution

It presents a novel TVD-based numerical method with hybrid line splittings tailored for EHL problems, improving efficiency and robustness over existing schemes.

Findings

Reduces computational complexity to O(n log n)

Handles high-order discretizations directly

Maintains robustness across a wide load parameter range

Abstract

In this study, we propose a class of total variation diminishing (TVD) schemes for solving pseudo-monotone variational inequality arises in elasto-hydrodynamic lubrication point contact problem. A limiter based stable hybrid line splittings are introduced on hierarchical multi-level grid. These hybrid splittings are designed by use of diffusive coefficient and mesh dependent switching parameter in the computing domain of interest. The spectrum of illustrated splittings is derived with the help of well known local Fourier analysis (LFA). Numerical tests validate the performance of scheme and its competitiveness to the previous existing schemes. Advantages of proposed splittings are observed in the sense that it reduces computational complexity (up to ($O(n\log n)$) and solve high order discretization directly (no defect-correction tool require) without perturbing the robustness of the solution procedure (i.e. it works well for large range of load parameters).

Tables

Table 1. Table 1 : Comparison of L ∞ superscript 𝐿 L^{\infty} -, L 1 superscript 𝐿 1 L^{1} -, and L 2 superscript 𝐿 2 L^{2} -error obtained for splitting L ​ s ​ 1 𝐿 𝑠 1 Ls1 in case of the linear convection-diffusion equation (Example  1 ϵ = 10 − 6 , κ = 1 / 3 formulae-sequence italic-ϵ superscript 10 6 𝜅 1 3 \epsilon=10^{-6},\kappa=1/3 ) over the domain Ω = [ − 1 , 1 ] × [ − 1 , 1 ] Ω 1 1 1 1 \Omega=[-1,1]\times[-1,1] .

$N$	$L^{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
16$\times$16	1.19566e-02	–	2.25624e-03	–	1.83208e-02	–
32$\times$32	2.62647e-03	2.1866	3.57540e-04	2.6577	2.92872e-03	2.6451
64$\times$64	5.70763e-04	2.2022	4.33084e-05	3.0454	3.64904e-04	3.0047
128$\times$128	1.06927e-04	2.4163	5.45271e-06	2.9896	4.73857e-05	2.9450
256$\times$256	1.92096e-05	2.4767	6.79793e-07	3.0038	6.09179e-06	2.9595
512$\times$512	3.40453e-06	2.4963	8.44721e-08	3.0085	7.74616e-07	2.9753

Table 2. Table 2 : Comparison of L ∞ superscript 𝐿 L^{\infty} -, L 1 superscript 𝐿 1 L^{1} -, and L 2 superscript 𝐿 2 L^{2} -error obtained for splitting L ​ s ​ 1 𝐿 𝑠 1 Ls1 in case of the linear convection-diffusion equation (Example  1 ϵ = 10 − 6 , κ = 0.0 formulae-sequence italic-ϵ superscript 10 6 𝜅 0.0 \epsilon=10^{-6},\kappa=0.0 ) over the domain Ω = [ − 1 , 1 ] × [ − 1 , 1 ] Ω 1 1 1 1 \Omega=[-1,1]\times[-1,1] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.27672e-02	–	1.49677e-03	–	1.36680e-02	–
$32\times 32$	2.73792e-03	2.2213	1.80364e-04	3.0529	1.82037e-03	2.9085
$64\times 64$	6.22587e-04	2.1367	7.33006e-05	1.2990	6.63061e-04	1.4570
$128\times 128$	2.07084e-04	1.5881	2.37525e-05	1.6257	2.13343e-04	1.6360
$256\times 256$	5.98623e-05	1.7905	6.73718e-06	1.8179	5.99206e-05	1.8321
$512\times 512$	1.58405e-05	1.9180	1.79203e-06	1.9106	1.58361e-05	1.9198

Table 3. Table 3 : Comparison of L ∞ superscript 𝐿 L^{\infty} -, L 1 superscript 𝐿 1 L^{1} -, and L 2 superscript 𝐿 2 L^{2} -error obtained for splitting L ​ s ​ 1 𝐿 𝑠 1 Ls1 in case of the linear convection-diffusion equation (Example  1 ϵ = 10 − 6 , κ = − 1.0 formulae-sequence italic-ϵ superscript 10 6 𝜅 1.0 \epsilon=10^{-6},\kappa=-1.0 ) over the domain Ω = [ − 1 , 1 ] × [ − 1 , 1 ] Ω 1 1 1 1 \Omega=[-1,1]\times[-1,1] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	2.32470e-03	–	1.56030e-02	–	2.05497e-02	–
$32\times 32$	1.01308e-03	1.1983	1.00995e-02	0.62754	9.31777e-03	1.1411
$64\times 64$	3.78032e-04	1.4222	4.35094e-03	1.2149	3.42296e-03	1.4447
$128\times 128$	1.07979e-04	1.8078	1.44691e-03	1.5884	9.67504e-04	1.8229
$256\times 256$	2.86739e-05	1.9129	4.47319e-04	1.6936	2.54980e-04	1.9239
$512\times 512$	7.39007e-06	1.9561	1.28974e-04	1.7942	6.53620e-05	1.9639

Table 4. Table 4 : Comparison of L ∞ superscript 𝐿 L^{\infty} -, L 1 superscript 𝐿 1 L^{1} -, and L 2 superscript 𝐿 2 L^{2} -error obtained for splitting L ​ s ​ 0 𝐿 𝑠 0 Ls0 in case of the linear convection-diffusion equation (Example  1 ϵ = 10 − 6 , κ = − 1.0 formulae-sequence italic-ϵ superscript 10 6 𝜅 1.0 \epsilon=10^{-6},\kappa=-1.0 ) over the domain Ω = [ − 1 , 1 ] × [ − 1 , 1 ] Ω 1 1 1 1 \Omega=[-1,1]\times[-1,1] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.62417e-02	–	2.30732e-03	–	2.04151e-02	–
$32\times 32$	1.01696e-02	0.67544	1.03223e-03	1.1605	9.53668e-03	1.0981
$64\times 64$	3.89903e-03	1.3831	3.65800e-04	1.4966	3.32527e-03	1.5200
$128\times 128$	1.20459e-03	1.6946	1.05973e-04	1.7874	9.49264e-04	1.8086
$256\times 256$	3.42856e-04	1.8129	2.84576e-05	1.8968	2.52429e-04	1.9109
$512\times 512$	9.05700e-05	1.9205	7.36748e-06	1.9496	6.49791e-05	1.9578

Table 5. Table 5 : Comparison of L ∞ superscript 𝐿 L^{\infty} -, L 1 superscript 𝐿 1 L^{1} -, and L 2 superscript 𝐿 2 L^{2} -error obtained for splitting L ​ s ​ 0 𝐿 𝑠 0 Ls0 in case of the linear convection-diffusion equation (Example  1 ϵ = 10 − 6 , κ = 1 / 3 formulae-sequence italic-ϵ superscript 10 6 𝜅 1 3 \epsilon=10^{-6},\kappa=1/3 ) over the domain Ω = [ − 1 , 1 ] × [ − 1 , 1 ] Ω 1 1 1 1 \Omega=[-1,1]\times[-1,1] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.17579e-02	–	2.25826e-03	–	1.83078e-02	–
$32\times 32$	1.76038e-03	2.7397	3.32640e-04	2.7632	2.72048e-03	2.7505
$64\times 64$	2.57573e-04	2.7728	4.20422e-05	2.9841	3.43166e-04	2.9869
$128\times 128$	3.47087e-05	2.8916	5.37451e-06	2.9676	4.37263e-05	2.9723
$256\times 256$	4.54820e-06	2.9319	6.78313e-07	2.9861	5.51381e-06	2.9874
$512\times 512$	6.02630e-07	2.9160	8.51091e-08	2.9946	6.91644e-07	2.9949

Table 6. Table 6 : Comparison of L ∞ superscript 𝐿 L^{\infty} -, L 1 superscript 𝐿 1 L^{1} -, and L 2 superscript 𝐿 2 L^{2} -error obtained for splitting L ​ s ​ 0 𝐿 𝑠 0 Ls0 in case of the linear convection-diffusion equation (Example  1 ϵ = 10 − 6 , κ = 0.0 formulae-sequence italic-ϵ superscript 10 6 𝜅 0.0 \epsilon=10^{-6},\kappa=0.0 ) over the domain Ω = [ − 1 , 1 ] × [ − 1 , 1 ] Ω 1 1 1 1 \Omega=[-1,1]\times[-1,1] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.24738e-02	–	1.50246e-03	–	1.36309e-02	–
$32\times 32$	2.00172e-03	2.6396	1.73122e-04	3.1175	1.68694e-03	3.0144
$64\times 64$	5.88728e-04	1.7656	6.97083e-05	1.3124	6.33640e-04	1.4127
$128\times 128$	1.84579e-04	1.6734	2.31011e-05	1.5934	2.08905e-04	1.6008
$256\times 256$	5.06126e-05	1.8667	6.64633e-06	1.7973	5.93352e-05	1.8159
$512\times 512$	1.28329e-05	1.9796	1.78059e-06	1.9002	1.57608e-05	1.9125

Table 7. Table 7 : Minimum film thickness result ( M = 20 , L = 10 formulae-sequence 𝑀 20 𝐿 10 M=20,L=10 ) for defect-correction κ = 0.0 𝜅 0.0 \kappa=0.0

Level	$H_m$	$H_m$(Moes)	$H_c$	$H_c$(Moes)	$H_c$(Moes)($p_x=0$)
1	1.99302e-01	1.92424	2.98940e-01	2.88624	2.77154
2	2.59716e-01	2.50753	3.70695e-01	3.57903	3.57760
3	2.70939e-01	2.61589	3.89566e-01	3.76122	3.75880
4	2.74629e-01	2.65151	3.94288e-01	3.80681	3.80443
5	2.75320e-01	2.65819	3.95428e-01	3.81782	3.81582
6	2.75525e-01	2.66016	3.95886e-01	3.82224	3.82034
7	2.75586e-01	2.66075	3.95962e-01	3.82297	3.82117

Table 8. Table 8 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} relative errors obtained with κ = 0.0 𝜅 0.0 \kappa=0.0 by Defect-Correction over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.57629e-01	–	4.56501e-03	–	9.85013e-02	–
$32\times 32$	1.75975e-01	-0.15884	2.01928e-03	1.1768	5.98804e-02	0.71806
$64\times 64$	1.69726e-01	0.052163	9.26960e-04	1.1233	3.78143e-02	0.66315
$128\times 128$	1.18555e-01	0.51765	3.56082e-04	1.3803	1.79500e-02	1.0749
$256\times 256$	7.20097e-02	0.71929	1.26752e-04	1.4902	7.87096e-03	1.1894
$512\times 512$	3.16527e-02	1.1859	4.43601e-05	1.5147	2.76403e-03	1.5098

Table 9. Table 9 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} relative errors obtained with κ = 1 / 3 𝜅 1 3 \kappa=1/3 by Defect-Correction over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.57629e-01	–	4.56501e-03	–	9.85013e-02	–
$32\times 32$	1.75975e-01	-0.15884	2.01928e-03	1.1768	5.98804e-02	0.71806
$64\times 64$	1.69726e-01	0.052163	9.26960e-04	1.1233	3.78143e-02	0.66315
$128\times 128$	1.18555e-01	0.51765	3.56082e-04	1.3803	1.79500e-02	1.0749
$256\times 256$	7.20097e-02	0.71929	1.26752e-04	1.4902	7.87096e-03	1.1894
$512\times 512$	3.16527e-02	1.1859	4.43601e-05	1.5147	2.76403e-03	1.5098

Table 10. Table 10 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained with κ = − 1.0 𝜅 1.0 \kappa=-1.0 by Defect-Correction over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.57629e-01	–	4.56501e-03	–	9.85013e-02	–
$32\times 32$	1.75975e-01	-0.15884	2.01928e-03	1.1768	5.98804e-02	0.71806
$64\times 64$	1.69726e-01	0.052163	9.26960e-04	1.1233	3.78143e-02	0.66315
$128\times 128$	1.18555e-01	0.51765	3.56082e-04	1.3803	1.79500e-02	1.0749
$256\times 256$	7.20097e-02	0.71929	1.26752e-04	1.4902	7.87096e-03	1.1894
$512\times 512$	3.16527e-02	1.1859	4.43601e-05	1.5147	2.76403e-03	1.5098

Table 11. Table 11 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained (M=20,L=10 case) with κ = 0.0 𝜅 0.0 \kappa=0.0 by splitting L h ​ s ​ 1 subscript 𝐿 ℎ 𝑠 1 L_{hs1} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$32\times 32$	7.99935e-02	–	3.25500e-03	–	4.31253e-02	–
$64\times 64$	6.76884e-02	0.240974	4.20806e-04	2.951430	1.35161e-02	1.673856
$128\times 128$	3.53135e-02	0.938689	1.14226e-04	1.881264	5.18955e-03	1.380998
$256\times 256$	1.01542e-02	1.798143	3.02821e-05	1.915354	1.35755e-03	1.934604
$512\times 512$	1.98897e-03	2.351983	8.51309e-06	1.830711	3.06834e-04	2.145475
$1024\times 1024$	4.02685e-04	2.304298	3.13898e-06	1.439387	8.16286e-05	1.910312

Table 12. Table 12 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained for (M=20,L=10 case) with κ = 1 / 3 𝜅 1 3 \kappa=1/3 by splitting L h ​ s ​ 1 subscript 𝐿 ℎ 𝑠 1 L_{hs1} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$32\times 32$	1.28495e-01	–	3.46499e-03	–	4.97302e-02	–
$64\times 64$	6.61681e-02	0.957504	4.17570e-04	3.052761	1.40651e-02	1.822002
$128\times 128$	3.34724e-02	0.983164	1.07470e-04	1.958084	5.05401e-03	1.476619
$256\times 256$	8.88278e-03	1.913889	2.70266e-05	1.991482	1.23452e-03	2.033478
$512\times 512$	1.64936e-03	2.429105	7.15546e-06	1.917264	2.47734e-04	2.317086
$1024\times 1024$	2.79280e-04	2.562122	2.77208e-06	1.368076	6.00344e-05	2.044930

Table 13. Table 13 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained (M=20,L=10 case) with κ = − 1.0 𝜅 1.0 \kappa=-1.0 by splitting L h ​ s ​ 1 subscript 𝐿 ℎ 𝑠 1 L_{hs1} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$32\times 32$	7.50604e-02	–	2.97122e-03	–	4.14394e-02	–
$64\times 64$	7.55099e-02	-0.008614	5.91844e-04	2.327767	1.69667e-02	1.288297
$128\times 128$	4.53322e-02	0.736130	1.91253e-04	1.629735	7.61954e-03	1.154930
$256\times 256$	1.61611e-02	1.488011	5.75179e-05	1.733400	2.50645e-03	1.604059
$512\times 512$	4.50872e-03	1.841736	1.67111e-05	1.783204	6.94586e-04	1.851420
$1024\times 1024$	1.10782e-03	2.024994	5.21125e-06	1.681105	1.89643e-04	1.872867

Table 14. Table 14 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained for (M=20, L=10) with κ = 0.0 𝜅 0.0 \kappa=0.0 by splitting L h ​ s ​ 2 subscript 𝐿 ℎ 𝑠 2 L_{hs2} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$32\times 32$	7.91753e-02	–	3.24093e-03	–	4.29201e-02	–
$64\times 64$	6.76405e-02	0.227163	4.21527e-04	2.942711	1.35422e-02	1.664191
$128\times 128$	3.53098e-02	0.937819	1.14185e-04	1.884252	5.18823e-03	1.384148
$256\times 256$	1.01543e-02	1.797978	3.02794e-05	1.914965	1.35750e-03	1.934290
$512\times 512$	1.99380e-03	2.348498	8.51277e-06	1.830636	3.07193e-04	2.143735
$1024\times 1024$	4.04313e-04	2.301976	3.13219e-06	1.442457	8.15121e-05	1.914059

Table 15. Table 15 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained for (M=20, L=10) with κ = 1 / 3 𝜅 1 3 \kappa=1/3 by splitting L h ​ s ​ 2 subscript 𝐿 ℎ 𝑠 2 L_{hs2} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$32\times 32$	1.27894e-01	–	3.45271e-03	–	4.95561e-02	–
$64\times 64$	6.61606e-02	0.950904	4.17669e-04	3.047297	1.40784e-02	1.815579
$128\times 128$	3.34692e-02	0.983138	1.07437e-04	1.958869	5.05304e-03	1.478260
$256\times 256$	8.88371e-03	1.913600	2.70267e-05	1.991034	1.23467e-03	2.033026
$512\times 512$	1.65390e-03	2.425290	7.15902e-06	1.916551	2.48217e-04	2.314452
$1024\times 1024$	2.80907e-04	2.557708	2.76808e-06	1.370876	5.99858e-05	2.048909

Table 16. Table 16 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained for (M=20, L=10) with κ = − 1.0 𝜅 1.0 \kappa=-1.0 by splitting L h ​ s ​ 2 subscript 𝐿 ℎ 𝑠 2 L_{hs2} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$32\times 32$	7.47880e-02	–	2.95607e-03	–	4.12735e-02	–
$64\times 64$	7.54384e-02	-0.012492	5.94019e-04	2.315099	1.70337e-02	1.276824
$128\times 128$	4.53370e-02	0.734610	1.91320e-04	1.634521	7.62081e-03	1.160376
$256\times 256$	1.61613e-02	1.488146	5.75274e-05	1.733667	2.50667e-03	1.604172
$512\times 512$	4.51054e-03	1.841171	1.67195e-05	1.782718	6.94549e-04	1.851624
$1024\times 1024$	1.10616e-03	2.027740	5.21053e-06	1.682030	1.89516e-04	1.873757

Table 17. Table 17 : Comparison of L ∞ superscript 𝐿 L^{\infty} , L 1 superscript 𝐿 1 L^{1} and L 2 superscript 𝐿 2 L^{2} errors obtained for EHL M=50 and L=10 with κ = 0.0 𝜅 0.0 \kappa=0.0 by splitting L h ​ s ​ 2 subscript 𝐿 ℎ 𝑠 2 L_{hs2} over the domain Ω = [ − 2.5 , 2.5 ] × [ − 2.5 , 2.5 ] Ω 2.5 2.5 2.5 2.5 \Omega=[-2.5,2.5]\times[-2.5,2.5] .

$N$	$L_{\infty }$-error	$p_{\infty }$	$L^1$-error	$p_1$	$L^2$-error	$p_2$
$16\times 16$	1.58602e-01	–	1.03810e-02	–	1.33934e-01	–
$32\times 32$	1.37546e-01	0.205497	2.42128e-03	2.100104	6.26015e-02	1.097253
$64\times 64$	9.91830e-02	0.471749	1.00043e-03	1.275150	3.00041e-02	1.061038
$128\times 128$	1.28502e-01	-0.373626	5.50322e-04	0.862272	2.15520e-02	0.477338
$256\times 256$	8.01042e-02	0.681841	3.32311e-04	0.727742	1.01793e-02	1.082183
$512\times 512$	4.33380e-02	0.886245	9.52456e-05	1.802810	3.69633e-03	1.461473

Equations

$$TV(u^n)=\sum _{-\infty }^{\infty }|u_{j+1}^n-u_j^n|=\sum _{-\infty }^{\infty }|{\Delta }_{j+1/2}{u}^n|$$

$${\Delta }_{j+1/2}{u}^n=u_{j+1}^n-u_j^n$$

$$u_t+f(u{)}_x=0.$$

$$TV(u^{n+1})\le TV(u^n)$$

$$u_i^{n+1}=u_i^n-c_i^n(u_i^n-u_{i-1}^n)+d_i^n(u_{i+1}^n-u_i^n)$$

$$c_i^n\ge 0,d_i^n\ge 0,c_i^n+d_i^n\le 1\forall i$$

$$\displaystyle\frac{\partial}{\partial x}\Big{(}\epsilon\frac{\partial u}{\partial x}\Big{)}+\frac{\partial}{\partial y}\Big{(}\epsilon\frac{\partial u}{\partial y}\Big{)}\leq\frac{\partial(\rho\mathscr{H})}{\partial x}\quad\in\quad\Omega$$

$$u\ge 0\in \Omega$$

$$\displaystyle u.\Big{[}\frac{\partial}{\partial x}\Big{(}\epsilon\frac{\partial u}{\partial x}\Big{)}+\frac{\partial}{\partial y}\Big{(}\epsilon\frac{\partial u}{\partial y}\Big{)}-\frac{\partial(\rho\mathscr{H})}{\partial x}\Big{]}=0\quad\in\quad\Omega,$$

$$u=0\text{on}\partial \Omega .$$

$$ϵ=\frac{\rho {\mathcal{H}}^3}{\eta \lambda },$$

$$\lambda =\frac{6{\eta }_0{u}_s{R}^2}{a^3{p}_H}.$$

$$\displaystyle\eta(u)=\exp\Bigg{\{}\Bigg{(}\dfrac{\alpha p_{0}}{z}\Bigg{)}\Bigg{(}-1+\Big{(}1+\dfrac{{u}p_{H}}{p_{0}}\Big{)}^{z}\Bigg{)}\Bigg{\}}.$$

$$\rho (u)=\frac{0.59\times 10^9+1.34u{p}_H}{0.59\times 10^9+u{p}_H}.$$

$$\mathcal{H}(x,y)={\mathcal{H}}_{00}+\frac{x^2}{2}+\frac{y^2}{2}+\frac{2}{{\pi }^2}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\frac{u(x^{{}^{\prime }},y^{{}^{\prime }})d{x}^{{}^{\prime }}d{y}^{{}^{\prime }}}{\sqrt{(x-x^{{}^{\prime }}{)}^2+(y-y^{{}^{\prime }}{)}^2}},$$

$${\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }u(x^{\prime },y^{\prime })d{x}^{\prime }d{y}^{\prime }=\frac{3\pi }{2}$$

$${\mathcal{D}}_f(x,y)=\frac{2}{{\pi }^2}\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}\frac{u(x^{{}^{\prime }},y^{{}^{\prime }})}{\sqrt{(x-x^{{}^{\prime }}{)}^2+(y-y^{{}^{\prime }}{)}^2}}d{x}^{{}^{\prime }}d{y}^{{}^{\prime }}.$$

$$\displaystyle\Omega^{hh}=\Big{\{}(x,y)\in\mathbb{R}^{2}\Big{|}x_{i^{{}^{\prime}}}-\frac{h}{2}\leq x\leq x_{i^{{}^{\prime}}}+\frac{h}{2},y_{j^{{}^{\prime}}}-\frac{h}{2}\leq y\leq y_{j^{{}^{\prime}}}+\frac{h}{2}\Big{\}}.$$

$${{\mathcal{D}}_f}_{i,j}={\mathcal{D}}_f(x_i,y_j)\approx \frac{2}{{\pi }^2}\sum _{i^{\prime }=0}^{n_x}\sum _{j^{\prime }=0}^{n_y}{{\mathcal{G}}^{hh}}_{i,i^{{}^{\prime }},j,j^{{}^{\prime }}}{u}_{i^{\prime },j^{\prime }}^{hh},$$

$${\mathcal{G}}_{i,i^{{}^{\prime }},j,j^{{}^{\prime }}}^{hh}=\underset{x_{i^{{}^{\prime }}}-\frac{h}{2}}{\overset{x_{i^{{}^{\prime }}}+\frac{h}{2}}{\int }}\underset{y_{j^{{}^{\prime }}}-\frac{h}{2}}{\overset{y_{j^{{}^{\prime }}}+\frac{h}{2}}{\int }}\frac{1}{\sqrt{(x-x^{{}^{\prime }}{)}^2+(y-y^{{}^{\prime }}{)}^2}}d{x}^{{}^{\prime }}d{y}^{{}^{\prime }}$$

$$x_i<x_{i^{{}^{\prime }}},x_i>x_{i^{{}^{\prime }}},x_i=x_{i^{{}^{\prime }}}\text{ and }{y}_j<y_{j^{{}^{\prime }}},y_j>y_{j^{{}^{\prime }}},y_j=y_{j^{{}^{\prime }}}$$

$$\displaystyle\mathscr{G}^{hh}_{i,i^{{}^{\prime}},j,j^{{}^{\prime}}}=\frac{2}{\pi^{2}}\Big{\{}|x_{+}|\sinh^{-1}(\frac{y_{+}}{x_{+}})+|y_{+}|\sinh^{-1}(\frac{x_{+}}{y_{+}})-|x_{-}|\sinh^{-1}(\frac{y_{+}}{x_{-}})$$

$$-|y_{+}|{\sinh }^{-1}(\frac{x_{-}}{y_{+}})-|x_{+}|{\sinh }^{-1}(\frac{y_{-}}{x_{+}})-|y_{-}|{\sinh }^{-1}(\frac{x_{+}}{y_{-}})$$

$$\displaystyle+|x_{-}|\sinh^{-1}(\frac{y_{-}}{x_{-}})+|y_{-}|\sinh^{-1}(\frac{x_{-}}{y_{-}})\Big{\}},$$

$$x_{+}=x_i-x_{i^{{}^{\prime }}}+\frac{h}{2},x_{-}=x_i-x_{i^{{}^{\prime }}}-\frac{h}{2}$$

$$y_{+}=y_j-y_{j^{{}^{\prime }}}+\frac{h}{2},y_{-}=y_j-y_{j^{{}^{\prime }}}-\frac{h}{2}.$$

$${\mathcal{H}}_{i,j}^{hh}:={\mathcal{H}}_{00}+\frac{x_i^2}{2}+\frac{y_j^2}{2}+\sum _{i^{\prime }}\sum _{j^{\prime }}{\mathcal{G}}_{|i-i^{\prime }|,|j-j^{\prime }|}^{hh}{u}_{i^{{}^{\prime }},j^{{}^{\prime }}}^{hh}=\text{H}{\mathcal{F}}_{i,j}^h,$$

$${\sigma }_{i,j}^h=\frac{r_{i,j}^h}{{\mathcal{G}}_{0,0}^{hh}},$$

$${r_{\text{J}}^h}_{i,j}=\text{H}{\mathcal{F}}_{i,j}^h-{\mathcal{H}}_{00}-\frac{x_i^2}{2}-\frac{y_j^2}{2}-\sum _{i^{\prime }}\sum _{j^{\prime }}{\mathcal{G}}_{|i-i^{\prime }|,|j-j^{\prime }|}^{hh}{\tilde{u}}_{i,j}^h$$

$${r_{GS}}_{i,j}^h=\text{H}{\mathcal{F}}_{i,j}^h-{\mathcal{H}}_{00}-\frac{x_i^2}{2}-\frac{y_j^2}{2}$$

$$-\sum _{i^{\prime }<i}\sum _{j^{\prime }}{\mathcal{G}}_{|i-i^{\prime }|,|j-j^{\prime }|}^{hh}{\bar{u}}_{i,j}^h-\sum _{i^{\prime }=i}\sum _{j^{\prime }<j}{\mathcal{G}}_{|i-i^{\prime }|,|j-j^{\prime }|}^{hh}{\tilde{u}}_{i,j}^h$$

$$-\sum _{i^{\prime }=i}\sum _{j^{\prime }>=j}{\mathcal{G}}_{|i-i^{\prime }|,|j-j^{\prime }|}^{hh}{\tilde{u}}_{i,j}^h-\sum _{i^{\prime }>i}\sum _{j^{\prime }}{\mathcal{G}}_{|i-i^{\prime }|,|j-j^{\prime }|}^{hh}{\tilde{u}}_{i,j}^h,$$

$${\tilde{\mathcal{G}}}_{i,j}^{hh}\simeq [{\mathcal{I}}_H^h{\mathcal{G}}_{i,.}^{hH}{]}_j,$$

$${\mathcal{W}}_i^h\simeq {\tilde{\mathcal{W}}}_i^h\stackrel{\text{def}}{=}{h}^d\sum _j{\tilde{\mathcal{G}}}_{i,j}^{hh}{u^{\ast }}_j^h=h^d\sum _j[{\mathcal{I}}_H^h{\mathcal{G}}_{i,.}^{hH}{]}_j{u^{\ast }}_j^h$$

Full text

1 Introduction

Elasto-hydrodynamic lubrication (EHL) is more often understood as a phenomenon of fluid film lubrication in which the natural process of hydrodynamic fluid film creation is governed due to deformation of contacting bodies and lubricant viscosity increases due to high pressure. Significant contributions have been made by many researchers in the development of more efficient and accurate methods for the study of EHL in last few decades (e.g.[17, 24, 10, 8, 6, 16, 1, 21, 18, 25]). It is well known that many numerical solutions of EHL model suffer lack of numerical stability and convergence during computation, if not tackled correctly. On the other hand, any direct solver such as Newton-Raphson technique takes a lot of computational storage and time (up to $O({n^{3}})$) to solve the model and hence it has no commercial use in practice. In 1992, Venner [24] has introduced a low order discretization for EHL model (see  1.1) which is stable for larger range of load parameters. However, author’s best knowledge stable schemes for the EHL model 1.1 are largely unavailable in literature which work well for very large range of load parameters other than Venner approach and in that sense it turns out to be a challenging problem in scientific community. The main numerical difficulty in these problems occurs due to lack of stable smoother and poor approximation of pressure profile near its steep gradient location by any standard iterative procedure. Also when applied load in contacting bodies are sufficiently high then many people observed wiggles in pressure and film thickness profile by using central or any high order scheme in convection term of Reynolds equation. One possible way to overcome the difficulty, people have used lower order discretization in convection term. In addition, for obtaining the high order stable, accurate solutions for such problems, researchers have applied lower order scheme in a defect corrected way through a suitable higher order discretization. However, such defect-correction [14, 13] setting most the time is not able to solve the difficulty in the sense that it does not reduce residual accurately due to poor conditioning of matrix in outer iteration (e.g.[20]). Furthermore, lower order schemes are more diffusive and allow to produce smoothing effect in the steep gradient region of solution and less accurate in the smooth part of the solution. This is the main motivation for present study to adopt total variation diminishing (TVD) approach for the EHL model problem. The reason behind TVD schemes for EHL model have been rarely applied so far in literature due to the fact that implementation is not obvious and straight forward as the case of linear-convection diffusion due to strong coupling of pressure and film thickness term in existing model. Therefore, in this article an attempt has been made to solve the problem generalizing TVD concept efficiently in the existing EHL model.

TVD schemes are understood as a generalized form of upwind based discretized schemes (more detailed definition will define later). Mostly, such schemes have been extensively devised for solving time dependent gas dynamics problems. Later on people have started to apply such concept for steady state problem in many CFD applications. Initially, the concept of TVD has been established by Harten and later by Sweby [11, 12, 22] to avoid unphysical wiggles in a numerical scheme. Harten also has given necessary and sufficient condition for a scheme to be TVD. To understand the concept, we first define the notation total variation $TV$ of a mesh function $u^{n}$ as

[TABLE]

having the following convention

[TABLE]

for any mesh function $u$ is used. Harten’s theory is understood in the form of conservation laws

[TABLE]

The numerical approximation of Eq. (3) is said to be TVD if

[TABLE]

Then Harten’s condition for any scheme to be TVD is explained below.

Theorem 1**.**

Let a general numerical scheme for conservation laws Eq. (3) is of the form

[TABLE]

over one time step, where the coefficients $c_{i}^{n}$ and $d_{i}^{n}$ are arbitrary value (In practice it may depend on values $u^{n}_{i}$ in some way i.e., the method may be nonlinear). Then $TV(u^{n+1})\leq TV(u^{n})$ provided the following conditions are satisfied

[TABLE]

There has been a very well developed TVD theory available in literature for time dependent problem. Additionally, this concept is also extended for steady state convection-diffusion case in the form of $M$- matrix [23] using appropriate flux limiting schemes [20, 14, 13, 19]. However, very little attention have been paid in developing TVD schemes for EHL problems. In this article, our aim to investigate a class of splitting for EHL model which is robust and high order accurate ( at least second order in smooth part of the solution ) for larger range of load parameters.

1.1 Model Problem

The following two dimensional circular point contact model problem is taken for numerical study defined below in the form of variational inequality written in non dimensional form

[TABLE]

where $u$ is non-dimensional pressure of liquid (lubricant) and $\Omega$ is sufficiently large bounded domain such that

[TABLE]

Here term $\epsilon$ is defined as

[TABLE]

where $\rho$ is dimensionless density of lubrication, $\eta$ is dimensionless viscosity of lubrication and speed parameter

[TABLE]

The non-dimensionless viscosity $\eta$ is defined according to

[TABLE]

Dimensionless density $\rho$ is given by

[TABLE]

The term film thickness $\mathscr{H}$ of lubricant is written as follows

[TABLE]

where $\mathscr{H}_{00}$ is an integration constant.

The dimensionless force balance equation is defined as follows

[TABLE]

All notations used in EHL model are defined in  A.

A schematic diagram of EHL point contact model is given in Fig. 1. Rest of the article is organized as followed. In Section. 2, few preliminaries are discussed which require in numerical study of EHL model which help in subsequent numerical analysis of the model. In Section 3, a series of splitting are constructed by imitating linear convection-diffusion model and linear EHL model. In Section 4, a hybrid splitting are constructed for solving our existing EHL model defined in Section 1. In Section 5, local Fourier analysis is performed to calculate quantitative estimate of splitting calculated in Section 3. In Section 6, numerical experiments are conducted to check the performance of present splitting and its improvement to EHL model. At the end of Section 7, overall conclusion is summarized.

2 Preliminaries

In this section, our main goal is to introduce few prerequisite theory which already used in our computation and cannot be ignored or avoided in the present analysis. Above nonlinear variational inequalities is solved numerically by using fixed point iteration theory [17, 3, 16]. The main challenge appears here in the form of producing a stable iterative smoother for EHL inequalities when the applied load on contacting bodies in EHL model become sufficiently large and after few iterations solution start blowing up. In such cases, iterative smoother for solving such model is stable only if nonlocal effect produced by film thickness equation is controlled by small change calculation in the iteration to make the overall effect local in updated pressure value. This effect is reduced by introducing special iterative smoother known as distributive smoother [24, 5, 4, 26]. The advantage of adopting such relaxation diminishes aggregation in film thickness computation and eventually leads to stable relaxation. Therefore, we need an extra care for computing film thickness term during each iteration. Let us define deformation integral $\mathscr{D}_{f}$ as

[TABLE]

We approximate the above integral Eqn. 14 taking pressure $u$ as piecewise constant function namely $u^{hh}_{i^{\prime},j^{\prime}}$ on sub-domain

[TABLE]

and discrete deformation

[TABLE]

where the coefficients $\mathscr{G}^{hh}_{i,i^{{}^{\prime}},j,j^{{}^{\prime}}}$ is written as

[TABLE]

and evaluated analytically. Above integration Eqn. 17 yields nine different results for the cases that are defined as

[TABLE]

respectively. The nine results are combined into one expression

[TABLE]

where

[TABLE]

Therefore film thickness in discretized form is written as

[TABLE]

where $\text{\tiny H}\mathscr{F}^{h}$ is right hand of the film thickness. For computing above discrete film thickness Eqn. 19, small change using relaxation is measured as

[TABLE]

where $\mathscr{G}^{hh}_{0,0}=\mathscr{G}^{hh}_{i=i^{\prime},j=j^{\prime}}$ and the residual ${r_{\text{J}}^{h}}_{i,j}$ for Jacobi relaxation is given by

[TABLE]

For Gauss-Seidel relaxation, residual ${r_{GS}}^{h}_{i,j}$ is given by

[TABLE]

where $\tilde{u}_{i,j}$ and $\bar{u}_{i,j}$ old and new updated values of pressure respectively.

2.0.1 Smooth kernel computation using MLMI

Suppose we want to solve integral of type Eqn. 19. If kernel $\mathscr{G}(x,y)$ is sufficiently smooth with respect to the variable $y$, we approximate discrete kernel $\mathscr{G}^{hh}_{i,j}$ by high order interpolation operator as

[TABLE]

where the high order interpolation operator is denoted by $\mathcal{I}^{h}_{H}$ and $\mathscr{G}^{hH}_{i,.}$ is injected from $\mathscr{G}^{hh}_{i,.}$ i.e., $\mathscr{G}^{hH}_{i,J}\stackrel{{\scriptstyle\text{def}}}{{=}}\mathscr{G}^{hh}_{i,2J}$. Superscript $h$ and $H$ denote the finer and the coarser grid respectively. Then the finer grid integral computation of Eqn. 19 is approximated on coarser grid in following way

[TABLE]

where

[TABLE]

Whenever kernel $\mathscr{G}(x,y)$ is also smooth enough with respect to $x$ variable, the discrete sum $\mathcal{W}^{h}_{i}$ is evaluated on coarse grid points $i=2I$ by use of high order interpolation operator $\mathcal{\hat{I}}^{h}_{H}$. It is written as

[TABLE]

where

[TABLE]

and where $\mathscr{G}^{HH}_{.,J}$ is injected from $\mathscr{G}^{hH}_{.,J}$, i.e., $\mathscr{G}^{HH}_{I,J}\stackrel{{\scriptstyle\text{def}}}{{=}}\mathscr{G}^{hH}_{2I,J}=\mathscr{G}^{hh}_{2I,2J}$.

2.0.2 Singular-Smooth or mild singular Kernel computation using MLMI

In general, kernel $\mathscr{G}$ has a mild singularity near a point $x=y$. We rewrite our coarse grid approximation by adding correction term near singularity in the following way (see [5])

[TABLE]

Since $\tilde{\mathscr{G}}^{hh}_{i,j}$ is an interpolation of $\mathscr{G}^{hh}_{i,j}$ itself using coarse grid points, the operator $(\mathscr{G}^{hh}_{i,j}-\tilde{\mathscr{G}}^{hh}_{i,j})$ is given by

[TABLE]

where $2p$ is the interpolation order and $\mathscr{G}^{2p}(\xi)$ is a $2p^{th}$ derivative of $\mathscr{G}$ at some intermediate point $\xi$. Thus if the derivative of $\mathscr{G}$ becomes small, the correction term become small and can be neglected. However, in case of singular smooth kernel ($i\simeq j$), we require the corrections in a neighborhood of $i=j(||j-i||\leq m\quad\text{or}\quad i-m\leq j\leq i+m)$. Thus Eq. (2.0.2) is simplified as follows

[TABLE]

Advantage of using multi-level procedure in film thickness $\mathscr{H}$ computation reduces integral complexity up to $O(n\log n)$. A schematic diagram of multi level multi integration procedure is given in Fig. 2.

2.1 Multi-Grid Method for variational inequality arising in EHL Problem

In this section, we discuss multi-grid method [9, 2] for variational inequality of EHL model. EHL problem is viewed as a linear complementarity problem [3, 19] of the form

[TABLE]

where $L$ is a linear differential operator. We want to solve the problem in discrete hierarchical sub-domains of the following form

[TABLE]

Hence discrete form of complementarity problem on level $l$ is written as

[TABLE]

Let $u_{l}$ and $v_{l}$ are an exact solution and approximated solution of above LCP Eqn. 33. Suppose that the error $e_{l}={u}_{l}-{v}_{l}$ is smooth after the iteration sweeping. Then complementarity problem satisfied for error equation $e_{l}$ on finer level is read as

[TABLE]

where residual $r_{l}=f_{1,l}-L_{l}v_{l}$. Such smooth error $e_{l}$ is approximated on a coarse grid without loosing any essential information. The LCP coarse grid equation for the coarse grid approximation of the error ${e}_{l-1}$ is therefore defined in PFAS by

[TABLE]

Since the problem is nonlinear and we are solving inequalities, we solve for full approximation ${v}_{l-1}={e}_{l-1}+I_{l}^{l-1}{v}_{l}$ but interpolate only $v_{l-1}$ back to fine grid. The main difference between multi-grid methods for equations and inequalities occur due to fact that, in case of fine grid converged solution $v_{l}=v_{l}^{*}$ the coarse grid correction equation should be zero. Consequently, we have the following relation

[TABLE]

(assume that operator $I^{l}_{l-1}$ keeps nonzero quantities nonzero).

Furthermore, for a converged solution of fine grid LCP problem the coarse grid correction provides us the following condition on restriction operators,

[TABLE]

Since $f_{1,l}-L_{l}v_{l}\equiv 0$ for any converge solution. Hence above inequalities 37 will satisfy for any rational choice of restriction operators $I^{l-1}_{l}$ and $\tilde{I}^{l-1}_{l}$. For capturing free boundary and for achieving fast convergence the bilinear interpolation operator $I_{l-1}^{l}$ is implemented only for unknowns on the inactive points that means,

[TABLE]

3 Linear study for convection-diffusion problem

Our specific interest in this Section is to develop an robust splitting for our EHL model. Such splitting is constructed by imitating series of linear model problem one by one. First we consider well known convection-diffusion problem of the form

Example 1**.**

[TABLE]

where $0<\epsilon<<1$ (note that we do not have any $y$ derivative in convection term). Then discretization of convective term for $(au)_{x}$ is performed as

[TABLE]

However, this scheme is only $O(h)$ accurate. Our interest here to increase accuracy at least smooth part without contaminating any wiggle in solution. Consider the Van Leer’s $\kappa$-schemes [15] for discretization term $(au)_{x}$ (for $a=\text{const}>0$) as

[TABLE]

(similar scheme can be constructed for $a<0$). The resulting discrete model Example. 1 by $\kappa$-scheme (take $\kappa=0$ here) is denoted by

[TABLE]

In general, above discrete equation. 1 do not produces $M$-matrix and many iterative splitting on $L_{\kappa}$ diverge. Therefore, this problem is solved using TVD scheme with help of appropriate flux limiters to prevent a solution from unwanted oscillation. Now consider $\kappa=-1$ then the second-order upwind scheme looks like ($a>0$)

[TABLE]

We enforce Eqn. 43 to satisfy TVD condition by multiply limiter functions in the additional terms $L_{\alpha},L_{\gamma}$ and $L_{\delta}$. Then following two type of discretization for convection term are presented here as

[TABLE]

and

[TABLE]

where $r_{i-1/2}=\dfrac{(u_{i+1,j}-u_{i,j})}{(u_{i,j}-u_{i-1,j})}$ and $r_{i-3/2}=\dfrac{(u_{i,j}-u_{i-1,j})}{(u_{i-1,j}-u_{i-2,j})}$.

In Fig. 3 represents graph of limiter function $(r,\phi(r))$ on which the resulting convection discretization term defined in Eqn. 43 and Eqn. 44 enforce to be TVD and higher order accurate (see [20]).

The discrete representation of Example 1 using Van-leer $\kappa$-scheme is defined as

[TABLE]

Moreover, in stencil notation it is represented as

[TABLE]

Then the discrete matrix equation $L_{\kappa}u=f$ is solved efficiently by the use of multi-grid. The related splitting is constructed by taking the matrix operator defined in Eqn. 47. In particular case, the splitting in $x$-direction is scanned as forward (or backward direction depending on flow direction) lexicographical order and it is represented as $S_{\kappa}=S_{\kappa}^{x_{f}}$ (or $S_{\kappa}^{x_{b}}$). For matrix operator $L_{\kappa}$, the forward splitting $S_{\kappa}^{x_{f}}$ is defined as

[TABLE]

where

[TABLE]

and therefore overall splitting is

[TABLE]

Now for a fixed $x$-line ($m$-grid points in $x$-direction) $(i,j_{0})_{(1\leq i\leq m)}$, we have the following

[TABLE]

$L^{0}_{\kappa}$ corresponds the operator to the unknowns $u^{*}$ which are scanned simultaneously. $L^{-}_{\kappa}$ corresponds the operator to the old approximation $u^{n}$, and $L^{+}_{\kappa}$ operator having updated values of $u^{n+1}$. Now by applying under-relaxation constant $\omega$ in above equation we have

[TABLE]

therfore splitting equation can be rewritten in corresponding change, $\sigma^{n+1}=u^{n+1}-u^{n}$ form as

[TABLE]

Now we construct series of splitting for solving Eqn. 1 as below.

Splitting : $L_{s0}$** This splitting is constructed by taking upwind operator $L_{1}$ plus a “positive” part of the second-order operators $L_{\alpha}$ and $L_{\beta}$ from Eqn. 45 and part of diffusion operator from Eqn. 47.**

[TABLE]

Splitting : $Ls1$** This splitting is constructed taking upwind operator $L_{1}$ plus a “positive” part of the second-order operators $L_{\alpha}$ from Eqn. 44 and part of diffusion operator from Eqn. 47.**

[TABLE]

Splitting : ${Ls2}$** In this case splitting coefficients $\mathcal{C}_{**}^{\kappa}$ correspond only to the first-order upwind operator $L_{1}$ of a discretized Eqn. 44 plus diffusion operator.**

[TABLE]

Splitting : $Ls3$** The third splitting named as $\kappa$- distributive line relaxation is constructed by assuming a ghost variable $\sigma_{*}$ (with the same cardinality as $\sigma$) such that $\sigma=\mathcal{D}\sigma_{*}$, where matrix $\mathcal{D}$ comes due to distributive change of the relaxation.i.e. We construct line-wise distributive splitting as**

[TABLE]

This splitting is understood in the following way: First, discretize Example 1 by $\kappa$-scheme and get the equation of the form as

[TABLE]

Now in the above splitting equation put the value of $u^{n+1}$ from Eqn. 51 and apply distributive splitting in the form of right preconditioner defined below.

[TABLE]

where the updated change in pressure and residual equation are denoted as

[TABLE]

**respectively. In other way, line distributive splitting consists of following two steps; In first step it calculates new ghost value approximation change $\sigma^{n+1}_{*}$. Second step calculates new approximation change $\sigma^{n+1}$.

Now applying above splitting along the $x$-direction in Example 1, the diffusive term is computed as**

[TABLE]

and convection term is computed as

[TABLE]

Other part of convective term which comes from van-leer discretization do not contain any distributive term as above explained and kept in right hand side during relaxation and overall splitting is written as follows

[TABLE]

after solving above equation for $\sigma$ along $x$ line direction updated solution $u^{n+1}$ is evaluated as

[TABLE]

**However, above splitting $Ls3$ Eqn. 54 is not robust and very rarely use in practice.

We are now interested in showing convergence of LCP through the above presented splitting. Let us consider domain $\Omega\in\mathbb{R}^{2}$ with boundary $\partial\Omega$, and consider known functions $f$ and $g$. Then find $u$ in a weak sense such that these inequalities hold**

Example 2**.**

[TABLE]

Therefore, discrete version of above problem (finite difference or finite volume) is written in the matrix form

[TABLE]

where $L$ is a $M$-matrix of order $m\times m$, $u$ and $f$ are $m\times 1$-column vector. It is well known that solving above discrete problem is equivalent to solving quadratic minimization problem of the form

[TABLE]

subjected to the constraints

[TABLE]

Theorem 2**.**

Let $u^{n}$ and $f^{n}$ are $m\times 1$-column vectors achieved by splitting algorithm (),*

[TABLE]

then we have $u^{n}\rightarrow u$ and $f^{n}\rightarrow f$ such that $u$ and $f$ is a solution of LCP problem.

Proof.

For the proof of this theorem we refer to see Cryer [7]. ∎

The following error estimates are easily established for LCP problem for algorithm described above.

Lemma 3**.**

Let $u$ is the exact solution of LCP problem define in Eqn. 3, also let $u^{n+1}$ is approximate solution obtained by the splitting of the form

[TABLE]

Then following conditions hold

[TABLE]

Proof.

Since From LCP problem we get

[TABLE]

and

[TABLE]

where

[TABLE]

Now consider the following LCP

[TABLE]

[TABLE]

[TABLE]

Now multiply $u^{T}$ in Eqn. 3 and combing with equality term we get

[TABLE]

similar way we also get

[TABLE]

Now by adding above two equations we get

[TABLE]

This implies that the following conditions hold

[TABLE]

[TABLE]

[TABLE]

Now rest of the proof is followed from Lemma 2.2 mentioned in [3]. ∎

Now we illustrate splitting for incompressible EHL model (we take $\rho,\eta$ and $\epsilon$ as constants here) in the form of inequalities as

Example 3**.**

[TABLE]

For incompressible EHL problem $\kappa$-line distributive Jacobi splitting is written as consider the convection term of above Example 3 as

[TABLE]

Now we will consider the following Splitting : ${Ls4}$

[TABLE]

Another possibility is to consider the following splitting as

[TABLE]

Hence overall equation is rewritten as Splitting : ${Ls5}$

[TABLE]

More general discussion on convergence of these splittings are given in Section 5.

4 TVD Implementation in Point Contact Model Problem

In this Section, we implement the splitting discussed in the last Section 3 and allow to extend it in EHL model. A hybrid splitting presented here and it is determined by measuring the value of

[TABLE]

This value is treated as switching parameter to perform two different splitting together while moving $x$ direction during the iteration. If the value

[TABLE]

then we apply line Gauss-Seidel splitting otherwise line Jacobi distributed splitting is incorporated in other words

[TABLE]

[TABLE]

These constructions are well justified as the region where $\epsilon$ tends to zero, we end up having an ill-conditioned matrix system in the form of dense kernel matrix appear in film thickness term. Therefore, distributive Jacobi line splitting is implemented as a right pre-conditioner to reduce the ill-conditioning of the matrix. However, in other part where $\epsilon$ is sufficiently large diffusion term dominates therefore we use Gauss line splitting. Considering the above setting in computational domain is quite demanding in EHL model as it allows us in reducing computational cost and storage issue. We replace $\kappa$ value in splitting constructed in Section 3 by incorporating appropriate limiter function $\phi$ there. In next section, we define these two splitting in more general form having limiter function involve in the splitting.

4.0.1 Limiter based Line Gauss-Seidel splitting

EHL point contact problem is solved in the form of LCP and therefore in this Section we seek an efficient splitting for Reynolds equation iterate along $x$-line direction to obtain the pressure solution. Now by using Theorem 2 and Lemma 3 we prove the convergence of the EHL solution. This splitting is explained in the following way: First calculate updated pressure in $x$-line direction as $\bar{u}_{i,j}=\tilde{u}_{i,j}+\sigma_{i}$ keeping $j$ fix at a time for all $j$ in $y$-direction and then apply change $\sigma_{i}$ immediately to update the pressure $\tilde{u}$. The successive pressure change $\sigma_{i}$ along the $x$-direction can be calculated as below

[TABLE]

where terms read as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

In above equation for each $i$,

[TABLE]

It is observed that the magnitude of the kernel $\mathscr{G}_{i,k,j,j}$ in equation  68 diminishes rapidly as distance $|k-i|$ increase and therefore, we avoid unnecessary computation expense by allowing value of $k$ up to three terms. So updated value of film thickness is rewritten as

[TABLE]

Hence, Eqn. (4.0.1) is illustrated as

[TABLE]

where $R_{i,j,\phi}$ and $\mathcal{C}_{i\pm.,\phi}$ are residual and coefficients of matrix arising due to linearized form involving the limiter function. This setting leads to a band matrix formulation which is solved using Gaussian elimination with minimum computational work ($O(n)$).

4.0.2 Limiter based Line-Distributed Jacobi splitting

The understanding philosophy of line distributed Jacobi splitting is more physical than mathematical. When diffusive coefficient tends to zero, pressure becomes large enough and non local effect of film thickness dominates in the region. Therefore a small deflection in pressure change produces high error in updated film thickness eventually leads blow up the solution after few iterations. This numerical instability is overcome by interacting with the neighborhood points during iteration. During this process the computed change of pressure at one point of the line are shared to its neighbor cells. In other words, a given point of a line new pressure $\bar{u}_{i,j}$ is computed from the summation of the changes coming from neighboring points plus the old approximated pressure $\tilde{u}_{i,j}$

[TABLE]

In this case, changes are incorporated only at the end of a complete iteration sweep. Therefore, overall splitting is derived as below

[TABLE]

The following notion used in Eqn. 4.0.2 defined as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

In above equation, discretization of convection term defined same as Line Gauss-Seidel relaxation case. However, due to distributive change of the pressure, the updated value of film thickness is described as

[TABLE]

where

[TABLE]

After few manipulation of Eqn. 4.0.2, we get system of band matrix which is solved using Gaussian elimination approach.

The force balance equation is incorporated in our numerical calculation by updating the constant value $\mathscr{H}_{00}$. The updated value of $\mathscr{H}_{00}$ is performed according to

[TABLE]

where $c$ is a relaxation parameter having range between $0.01-0.1$.

5 Fourier Analysis

Performance and asymptotic estimate of above splitting is measured through the Fourier analysis by considering infinite grid

[TABLE]

and infinite grid function defined on ${\mathbb{G}^{f}}_{h}$ by the linear span of the Fourier components

[TABLE]

These basis functions $e^{i{\xi\theta}}\in\mathbb{T}^{h}$ are orthogonal with respect to the inner product

[TABLE]

where $u_{h},v_{h}\in\mathbb{T}^{h}$. Furthermore, we will define orthogonal space to identity function $\mathbb{I}\in\mathbb{T}^{h}$ as

[TABLE]

Moreover, discrete solution $u_{h}$ is described as Fourier transform $\hat{u}$ a linear combinations of the basis functions $e^{i\xi\theta}\in\mathbb{T}^{h}$

[TABLE]

The Fourier space $\mathbb{T}^{h}$ is illustrated as four-dimensional subspaces

[TABLE]

[TABLE]

We say discretized PDE of the form

[TABLE]

is solvable if $f_{h}\in\mathbb{T}^{h}_{\bot}$. Moreover, solution will be unique if $u_{h}\in\mathbb{T}^{h}_{\bot}$. Let relaxation method defined via operator splitting as

[TABLE]

where $\tilde{u}_{h}$ and $\bar{u}_{h}$ are old and updated approximation to the solution $u_{h}$. Now we are interested in constructing a splitting which reduce our computed error significantly. Such behavior is investigated by measuring error equation as

[TABLE]

where $\tilde{e}_{h}=u_{h}-\tilde{u}_{h}$, $\bar{e}_{h}=u_{h}-\bar{u}_{h}$ and $\mathcal{S}_{h}:=-(L^{+}_{h})^{-1}L^{-}_{h}$. Now apply Fourier transform in above equation for $\hat{L}^{+}_{h}(\theta)\neq 0$ we have following relation

[TABLE]

and smoothing factor notation as

[TABLE]

where $\hat{\mathcal{S}}_{h}(\theta):=-\hat{L}^{-}_{h}(\theta)/\hat{L}^{+}_{h}(\theta)$.

5.0.1 Fourier analysis of $\kappa$ splitting

Let $\tilde{u}^{h}_{i,j}$ current updated to the solution for given $j$ line we are solving equations. For given $j$ a new updated $\bar{u}^{h}_{i,j}$ for all $i$ of that line according to

[TABLE]

for $2\leq i\leq(n_{x}-1)$ and for given value $j$ such that $1\leq j\leq n_{y}-1$ holds. During Gauss-Seidel line relaxation, we will use previously computed new solution of line $j-1$ in our next new updated solution of line $j$. Hence error equation is written as

[TABLE]

and $\kappa$-smoothing factor is denoted as

[TABLE]

where $\alpha_{1}=\epsilon/h^{2}$ and $\beta=a/h$. Smoothing factor plot is given in Fig. 4

Two grid iteration matrix is written as

[TABLE]

and two grid error equation is defined as

[TABLE]

Here by multiplying $C_{h}^{2h}$ to the space $\mathbb{T}^{h}_{\theta}$, where ${\theta}\in\tilde{\varTheta}^{00}=\varTheta^{00}-\{{\theta}:L_{2h}(2\theta^{00})=0\}$ leaves the space invariant.

[TABLE]

Fourier representation of two grid is performed in following way

[TABLE]

Spectral radius is computed in the following way

[TABLE]

where

[TABLE]

The Fourier symbols of the multi-grid operators for each harmonic in $\mathbb{T}^{h}_{\theta}$ is calculated as follows:

[TABLE]

For the transfer operators

[TABLE]

Since we can always get a nonsingular matrix $P$ same order as $C^{2h}_{h}$ such that $PC^{2h}_{h}P^{-1}=Q^{2h}_{h}$ holds, where $Q_{h}^{2h}$ a block matrix consisting of $4\times 4$ diagonal block $\tilde{Q}^{2h}_{h}(\theta)$ looks for all $\theta\in\tilde{\Theta}_{00}$ like

[TABLE]

then the smoothing factor is equivalent to

[TABLE]

Computation of $\mu$ is important for observing two-grid convergence during relaxation. In next Section we illustrate a criterion for two-grid convergence.

5.1 Convergence criterion of hybrid splitting

In this section, we give a general criteria for the convergence study of hybrid schemes used in our EHL model problem. Let us reconsider linear system

[TABLE]

where $[L_{\kappa}]_{m\times m}$ a regular matrix (for definition see ****[23]****) and $f$ and $u$ are known values. For applying hybrid splitting in above equation matrix $L_{\kappa}$ is understood as

[TABLE]

where $[L_{\kappa}^{\Omega_{\epsilon}}]$ and $[L_{\kappa}^{\Omega^{\prime}_{\epsilon}}]$ are regular applied splittings in

[TABLE]

and

[TABLE]

sub-domains respectively.

Now assume that $[L_{\kappa}^{\Omega_{\epsilon}}]$ has the following splitting

[TABLE]

where $M_{\kappa}^{\Omega_{\epsilon}}$ is a regular easily invertible matrix and $N_{\kappa}^{\Omega_{\epsilon}}$ is a positive rest matrix. Then our splitting can be defined as

[TABLE]

Then above iteration will converge for any initial guess $u^{0}$ if following theorem holds

Theorem 4**.**

Let $L_{\kappa}^{\Omega_{\epsilon}}=M_{\kappa}^{\Omega_{\epsilon}}-N_{\kappa}^{\Omega_{\epsilon}}$ be a regular splitting of matrix $L_{\kappa}^{\Omega_{\epsilon}}$ and $(L_{\kappa}^{\Omega_{\epsilon}})^{-1}\geq 0$, then we have

[TABLE]

Proof.

For the proof of this theorem we refer to see Varga [23]. ∎

Now we will prove other part of matrix splitting $L_{\kappa}^{\Omega^{\prime}_{\epsilon}}$. This part of matrix there is no straightforward splitting is available (see ****[23, 26]****). Let $L_{\kappa}^{\Omega^{\prime}_{\epsilon}}$ is regular, but dense and the designing suitable splitting in the sense of Varga is complicated. Suppose if it is possible to construct nonsingular matrix $L^{r}_{\kappa}$ such that equation below

[TABLE]

is easy to solve and we can rewrite splitting as

[TABLE]

Then for above splitting our iteration is denoted as

[TABLE]

Therefore above iteration will converge for any initial guess if following theorem holds

Theorem 5**.**

Let $(M_{\kappa}^{\Omega^{\prime}_{\epsilon}}-N_{\kappa}^{\Omega^{\prime}_{\epsilon}})(L^{r}_{\kappa})^{-1}$ be a regular splitting of matrix $L_{\kappa}^{\Omega^{\prime}_{\epsilon}}$ and $(L_{\kappa}^{\Omega^{\prime}_{\epsilon}})^{-1}\geq 0$, then we have

[TABLE]

The following theorem providing sufficient conditions for the convergence of the two-grid method $Q_{2}$ ( define in Eqn.105) is due to Hackbusch.

Theorem 6**.**

Let us assume that $\mathcal{S}_{l}$ is a smoothing operator for $K_{l}$ that means there exist $\eta(\nu)$ and $\nu^{\prime}(h)$ so that the following condition holds

[TABLE]

and also assume that operator $K_{l}$ is approximated accurately (by prolongation and restriction operator) in the following sense such that $\exists\quad C_{A}\rightarrow 0$, independent of $h$ so that

[TABLE]

then there exist $h$ and $\nu\in\mathbf{N}$:

[TABLE]

holds for $\nu$ with $\nu^{\prime}(h_{l})\geq\nu\geq\nu(h_{l})$ and $h_{2}\leq h$ and the two-grid method $Q_{2,l}$ from Eqn. 97 converges monotonically, independently of $h$.

Proof.

It follows straight way by taking $Q_{2,l}(\nu,0)=(K^{-1}_{l}-P(K_{l-1})^{-1}R)(K_{l}\mathcal{S}_{l}^{\nu})$. ∎

6 Numerical Results

In Section 3, we have illustrated TVD implementation for solving linear convection-diffusion problem through a class of splittings. Now we investigate the performance of mentioned splittings and compare the results with classical defect-correction. For numerical tests we consider analytical solution as $u=x^{4}+y^{4}$ from Oosterlee ****[20]****. All numerical computations is performed on author’s personal laptop having 2GB RAM and Intel(R) Core(TM) i3-2328M CPU @ 2.20GHz. Dirichlet boundary is imposed for all test cases on domain \Omega=\Big{\{}(x,y);-1\leq x\leq 1,-1\leq y\leq 1\Big{\}}. For all numerical experiments, we take diffusion coefficient $\epsilon=10^{-6}$ and $\kappa=-1.0,0.0,1/3$. Numerical tests are performed for the problem given as Example 1 using $Ls0$ splitting, $Ls1$ splitting and classical defect-correction technique using hierarchical multi-level grid. Computational results of relative error and corresponding order in $L^{1},L^{\infty},L^{2}$-norms are presented on Table 1- 6 on the finest grid level ($7^{th}$** level using $3V(2,1)$ cycle).**

$L^{2}$** norm error is evaluated in the following way**

[TABLE]

where $H$ is the mesh size on grid $k-1$, $\bar{u}^{k}$ is the converged solution on grid $k$ and $d$ denotes the dimension of the problem. The order of convergence is derived as

[TABLE]

where $p_{2}$ is the order of discretization in $L^{2}$ norm. We also calculate $L^{\infty}$ and $L^{1}$-error and corresponding order in similar fashion. From numerical experiments we observe that splitting $Ls0$ and $Ls1$ always show fast residual decay compare to classical defect-correction. Fig. 5 and Fig. 6 present the residual decay results for $Ls0$ splitting , $Ls1$ splitting and classical defect-correction technique for $\kappa=0.0,1/3$. Moreover, residual decay of splitting $Ls1$ is more better than splitting $Ls0$. On the other hand, we observe that splitting $Ls0$ has larger range of robustness ($-1.0\leq\kappa\leq 0.9$) than splitting $Ls1$ ($-1.0\leq\kappa\leq 0.8$).

6.1 Test case for numerical experiment of EHL problem

In this section, we perform numerical experiments on EHL model defined in Section 1. We take Moes (****[18]****) dimensionless parameters (which is denoted by $M$ and $L$), where $L$ is fixed at $10$ while $M$ is varied between $20-1000$. For all test cases, we fix the parameter $\alpha=1.7\times 10^{-8}$ over domain $\Omega=[-2.5,2.5]\times[-2.5,2.5]$. In all cases , we refine grid up to $(1024+1)\times(1024+1)$ points on finest level and coarse grid up to $(32+1)\times(32+1)$ points on the coarsest level (except extremely high load case we choose coarse grid ($64+1)\times(64+1)$). A class of limiter are applied to solve the problem discussed in Section 3 and 4. However, for checking performance of splittings, we use value $\kappa=0.0,1/3,-1.0$ in our numerical analysis. In Fig. 8, we represent film thickness profile $\mathscr{H}$ in inverted form. Four load cases (a)$M=20,L=10$, (b)$M=50,L=10$, (c) $M=100,L=10$ and (c) $M=1000,L=10$ are solved using the TVD schemes. The fully converged pressure as well as film thickness profiles and their plot results are represented in Fig. 8-Fig.12. Comparisons of relative error in $L^{2},L^{1}$ and $L^{\infty}$ norms between $\kappa$ splittings and defect correction schemes are performed which are presented in Table. 7- 17. Experimental results show that order of convergence of classical defect-correction is almost similar to splittings $L_{hs1}$ and $L_{hs2}$. However, splittings $L_{hs1}$ and $L_{hs2}$ have slightly better residual decay in comparison with classical defect-correction which can be seen in Fig. 7.

7 Conclusion

A limiter based hybrid line splittings have been outlined for solving EHL point contact problem (in the form of LCP) on hierarchical multi level grid. The key idea of using such splitting to facilitate artificial diffusion only the region of steep gradient of pressure profile and to improve the accuracy on the other part (smooth region of pressure profile) of the domain. These illustrated splittings have been devised by bringing left hand side matrix in $M$-matrix form using second order discretization of Reynolds equation and rest term on the right hand side. Additionally, the hybrid line splitting has been designed with help a switcher which depends upon magnitude of $\epsilon/h$. When $\epsilon/h\leq 0.6$, we have applied distributive Jacobi line splitting else, we have implemented Gauss-Seidel line splitting during updating new solution. The derived switcher is important as it noticeably allows us in reducing the ill-conditioning of the discretized matrix when $\epsilon$ is almost equal to zero. The robustness of the splittings have been analyzed performing series of numerical experiments. Moreover, robustness range of splittings has been investigated and compared with other splittings. For linear $\kappa-$ discretization, we have performed Fourier analysis in order to validate the multi-grid convergence behavior theoretically. Numerical experiments conform that the performance of these hybrid line splittings are robust not only for linear case but also for EHL model too. A remarkable achievement of these splittings are that it helps us in developing of higher-order discretization without losing stability in relaxation and without the use of double discretization scheme like defect-correction technique in multi-grid solver. Numerical experiments confirm that residual decay of direct splittings are comparably better than classical defect-correction. In this study, we have analyzed the performance of splittings through known limiters available in literature which works satisfactory in all study cases. Another remarkable advantage of the adopted splittings can be noted as it does not demand any extra tuning parameter and produces reasonable numerical solution for large range of load variation.

8 Acknowledgment

**This work is fully funded by DST-SERB Project reference no.PDF/2017/000202 under N-PDF fellowship program and working group at the Tata Institute of Fundamental Research, TIFR-CAM, Bangalore. Author is highly indebted to IIT Kanpur for all kind of support that facilitated the completion of this work. **

Appendix A Some Notation used in EHL model

$p_{H}\rightarrow$** Maximum Hertzian pressure.**

$\eta_{0}\rightarrow$** Ambient pressure viscosity.**

$H_{00}\rightarrow$** Central offset film thickness.**

$a\rightarrow$** Radius of point contact circle.**

$\alpha\rightarrow$** Pressure viscosity coefficient.**

$u_{s}=u_{1}+u_{2}$, where $u_{1}$ upper surface velocity and $u_{2}$ lower surface velocity respectively.

$p_{0}\rightarrow$** Constant ($p_{0}=1.98\times 10^{8}$), $z$ is pressure viscosity index ($z=0.68$).**

$R\rightarrow$** Reduced radius of curvature defined as $R^{-1}=R_{1}^{-1}+R_{2}^{-1}$,**

where $R_{1}$ and $R_{2}$ are curvature of upper contact surface and lower contact surface respectively.

$L$** and $M$ are Moes parameters and they are related as below.**

$L=G(2U)^{\frac{1}{4}},M=W(2U)^{-\frac{1}{2}}$**, where **

$2U=\dfrac{(\eta_{0}u_{s})}{(E^{{}^{\prime}}R)},W=\dfrac{F}{E^{\prime}R},p_{H}=\dfrac{(3F)}{(2\pi a^{2})}$.

$\sigma^{n+1}=u^{n+1}-u^{n}$** denote as difference between latest approximation solution $u^{n+1}$ and its predecessor $u^{n}$.**