"Icebreaker problem": Dividing set into partitions that maximize pairwise shared memberships?

Question

Suppose I have a set of N*L elements. I wish to generate a sequence of k partitions of size L such that I maximize the number of pairs of elements that share membership of the same partition. What kind of algorithm would optimize this?

To make it easier to understand, imagine an icebreaker situation: Suppose I have a new class of 15 students. Every week over the course of the first 1 month (4 weeks), I want to divide the students into 3 different teams of 5 so that by the end of the month, all of the students have worked with as many of their classmates as possible.

Naively choosing random subsets without replacement turns out to be a rather effective (and easily implementable) solution:

# -*- coding: utf-8 -*-

import random
import copy

import pandas as pd
import matplotlib.pyplot as plt

plt.style.use('ggplot')


NUM_TEAMS = 3                             # N
TEAM_SIZE = 5                             # L
CLASS_SIZE = NUM_TEAMS * TEAM_SIZE        # N*L
NUM_WEEKS = 4
NUM_PARTITIONS = NUM_TEAMS * NUM_WEEKS    # k

NUM_SIMULATIONS = 10000


def run_simulation(students, max_pairings):

    weeks = [[] for _ in range(NUM_WEEKS)]
    pairings = dict((student, set()) for student in students)

    for week_num, week in enumerate(weeks):

        unassigned_students = copy.deepcopy(students)
        teams = [[] for _ in range(NUM_TEAMS)]

        for idx, team in enumerate(teams):
            team = set(random.sample(unassigned_students, TEAM_SIZE))
            teams[idx] = team
            for student in team:
                pairings[student] = pairings[student].union(team)
            unassigned_students = list(set(unassigned_students) - team)

        weeks[week_num] = teams

    pairings = dict((k, len(pairings[k])) for k in pairings)
    total_pairings = sum(pairings[k] for k in pairings)

    if total_pairings > max_pairings:
        max_pairings = total_pairings
        for week in weeks:
            print(week)
        print("Pairings by member:", pairings)
        print("Total pairings: ", total_pairings, '\n')

    return max_pairings


if __name__ == '__main__':

    students = range(CLASS_SIZE)
    max_pairings = 0
    results = []

    for iteration in range(NUM_SIMULATIONS):
        max_pairings = run_simulation(students, max_pairings)
        results.append((iteration, max_pairings))

    df = pd.DataFrame(results, columns=['Iteration', 'Max pairings'])
    df.set_index('Iteration')
    
    df['Max pairings'].plot()
    plt.xlabel('Iteration')
    plt.ylabel('Max number of pairings')
    plt.show()

It's easy to extend this naive solution to a Monte Carlo simulation, which has rather fast convergence:

Solving the above parameterization, I can arrive at a schedule with 213 pairings within about 10^6 simulations, but it's very hard to come up with a solution that beats that...