I don't know whether there are any good heuristics, but I can suggest one approach: you could try to use an ILP solver on this task. Introduce zero-or-one variables $x_{i,v}$, $y_{i,u,v,w}$, where $x_{i,v}=1$ means that vertex $v$ is contained in partition $P_i$, and $y_{i,u,v,w}=1$ means that $u,v,w$ are a 3-clique in $P_i$. (Only introduce variables $y_{i,u,v,w}$ for triples of vertices $u,v,w$ that are a 3-clique in the original graph $G$; all others are effectively hardcoded at 0.) You can add constraints
$$y_{i,u,v,w} \ge x_{i,u} + x_{i,v} + x_{i,w} - 2,$$
$$y_{i,u,v,w} \le x_{i,u}, y_{i,u,v,w} \le x_{i,v}, y_{i,u,v,w} \le x_{i,w},$$
$$\sum_v x_{i,v} \le M, \sum_i x_{i,v}=1.$$
Then, maximize the objective function $\sum_{i,u,v,w} y_{i,u,v,w}$ subject to these constraints, using an off-the-shelf ILP solver.