Family of Generalized Integrals ${I}(a,b,p)=\int_0^{ab} \left( \left\{\frac{x}{a}\right\}-p\right) \left( \left\{\frac{x}{b}\right\}-p\right) \; dx$

This just addresses the simplest questions in the OP, namely the estimation of $I(a,b;p)$ when $d:=g.c.d(a,b)>1$. In such case, there are integers $q,r$ such that $a=qd$, $b=rd$ and $g.c.d(q,r)=1$

The change of variable $u=x/d$ and the periodicity of $x\mapsto\{x\}$ results in

$$ \begin{align} I(a,b;p)&=d\int^{dqr}_0\Big(\big\{\tfrac{x}{q}\big\}-p\Big)\Big(\big\{\tfrac{x}{r}\big\}-p\Big)\,dx=d\sum^{d-1}_{j=0}\int^{(j+1)qr}_{jqr}\Big(\big\{\tfrac{x}{q}\big\}-p\Big)\Big(\big\{\tfrac{x}{r}\big\}-p\Big)\,dx\\ &= d\sum^{d-1}_{j=0}\int^{qr}_0\Big(\big\{\tfrac{x+jqr}{q}\big\}-p\Big)\Big(\big\{\tfrac{x+jqr}{r}\big\}-p\Big)\,dx\\ &=d^2\int^{qr}_0\Big(\big\{\tfrac{x}{q}\big\}-p\Big)\Big(\big\{\tfrac{x}{r}\big\}-p\Big)\,dx \end{align}$$

Form the case of relative prime one obtained (the OP provides a sketch of the proof for this case):

$$\begin{align} \frac{1}{ab}I(a,b;p)= \frac{1}{4}{\left(1-2p\right)}^2+\frac{1}{12}\frac{g.c.d(a,b)}{l.c.m(a,b)} \end{align}$$

where $l.c.m(a,b)$ is the lowest common multiple of $a$ and $b$. When $p=1/2$ we recover Franel's formula.

For higher order integrals, I have a few references that may be useful to those interested in this question:

Franel integrals of order three

Franel Integrals of over four

Multiple Franel integrals

Here I provide further details for the expression $I(a,b;p)$ when $g.c.d(a, b)=1$.

Without loss of generality assume $a>b$. The $a=bq+r$ for $q\in\mathbb{N}$ and $1\leq r<b$

$$ \begin{align} I(a,b;p)&=\sum^{b-1}_{k=1}\int^{(k+1)a}_{ka}\Big(\big\{\tfrac{x}{a}\big\}-p\Big)\Big(\big\{\tfrac{x}{b}\big\}-p\Big)\,dx = \sum^{b-1}_{k=0}\int^a_0\Big(\big\{\tfrac{x+ak}{a}\big\}-p\Big)\Big(\big\{\tfrac{x+ak}{b}\big\}-p\Big)\,dx\\ &=\sum^{b-1}_{k=0}\int^a_0\Big(\big\{\tfrac{x}{a}\big\}-p\Big)\Big(\big\{\tfrac{x+rk}{b}\big\}-p\Big)\,dx=\sum^{b-1}_{k=0}\sum^{a-1}_{\ell=0}\int^{\ell+1}_{\ell}\Big(\big\{\tfrac{x}{a}\big\}-p\Big)\Big(\big\{\tfrac{x+rk}{b}\big\}-p\Big)\,dx\\ &=\sum^{b-1}_{k=0}\sum^{a-1}_{\ell=0}\int^1_0\Big(\big\{\tfrac{x+\ell}{a}\big\}-p\Big)\Big(\big\{\tfrac{x+\ell+rk}{b}\big\}-p\Big)\,dx\ \end{align} $$ Since $g.c.d(a,b)=g.c.d(b,r)=1$, for each $0\leq\ell<a$ fixed, $\ell+kr$ ranges over $\{0,\ldots,b-1\}\mod \,b$ as $k$ over $\{0,\ldots,b-1\}$. Hence $$ I(a, b;p)=\sum^{b-1}_{j=0}\sum^{a-1}_{\ell=0} \int^1_0\Big(\big\{\tfrac{x+\ell}{a}\big\}-p\Big)\Big(\big\{\tfrac{x+j}{b}\big\}-p\Big)\,dx=\sum^{b-1}_{j=0}\sum^{a-1}_{\ell=0}\int^1_0 \Big(\tfrac{x+\ell}{a}-p\Big)\Big(\tfrac{x+j}{b}-p\Big)\,dx $$

The rest is as the OP indicated.

I ignore thus far whether a similar argument carries over for higher orders under the assumption that $g.c.d(a_1,\ldots,a_n)=1$.

This answer is only for the case $\gcd(a_1, a_2), \gcd(a_1, a_3),... = 1$ (i.e. the $\gcd$ of any pair of $a_i$ is $1$).

Starting from what you have already done: $$I(a_1, a_2, ..., a_n, p) = \int_0^1 \prod_{i=1}^n \left(\frac{a_i-1}{2} +t-a_ip \right)dt$$

Let $c_i = \frac{a_i-1}{2}-a_ip$. Then the task is to find $$\int_0^1 \prod_{i=1}^n\left( t+c_i \right)dt$$

The integrand can be expanded as $$\int_0^1 \sum_{k=0}^nS_{k, n}x^{n-k} dt$$

where $S_{k, n}$ is the sum of the product of all $\binom{n}{k}$ "$k$-tuples" from $c_1, c_2, ..., c_n$ (except for $S_{0, n} = 1$). For example, $S_{2, 4} = c_1c_2+c_1c_3+c_1c_4 + c_2c_3+c_2c_4+c_3c_4$ and $S_{3, 4} = c_1c_2c_3+c_1c_2c_4+c_1c_3c_4+c_2c_3c_4$. Let $s_{k, n}$ be defined in a similar way, but for $a_i$ instead of $c_i$

Then the integral is $$\sum_{k=0}^{n} \frac{S_{k, n}}{n+1-k}$$

For $n = 3$, this gives $$\frac{(1-2p)^3a_1a_2a_3}{8}+\frac{(1-2p)(a_1+a_2+a_3)}{24} = $$ $$\frac{(1-2p)^3 s_{3, 3}}{8} + \frac{(1-2p) s_{1, 3}}{24}$$

For $n = 4$, this gives $$\frac{(1-2p)^4 s_{4, 4}}{16}+\frac{s_{2, 4}}{48}(1-2p)^2 + \frac{1}{80}$$

For $n = 5$, this gives $$\frac{(1-2p)^5 s_{5, 5}}{32} + \frac{(1-2p)^3 s_{3, 5}}{96} + \frac{(1-2p)s_{1, 5}}{160}$$

In general it seems like $$I(a_1, ..., a_n, p) = \sum_{1 \le k \le n+1, k\pmod2 = 1} \frac{(1-2p)^{n+1-k} s_{n+1-k, n}}{k\cdot 2^n}$$

although I have not confirmed this.