Для случая электростатики (неподвижные заряды и неизменные электрические поля) теорему Гаусса (физическую) можно доказать элементарными методами, то есть по методу Ньютона.
Рассмотрим каждый точечный заряд, создающий вокруг себя кулоновское электрическое поле

Суммарное электрическое поле - это поле, получаемое векторной суммой этих кулоновских полей, рассчитанной в каждой точке пространства (принцип суперпозиции). Соответственно, и интеграл от поля можно считать как сумму интегралов для каждого отдельного точечного заряда.
Дальше рассмотрим два случая:
I. Точечный заряд охвачен поверхностью интегрирования.
Разобьём всё пространство вокруг заряда на тонкие конические трубки прямыми лучами, которые исходят из точки заряда. Каждая такая трубка будет содержать телесный угол

а все вместе - они будут содержать телесный угол

Поверхность интегрирования будет пересекать каждую трубку 1 раз, а если она пересечёт какую-то трубку несколько раз, то это будет нечётное число раз, и ориентация пересечений будет разной:

раз нормаль к поверхности будет направлена "наружу", от заряда, и

раз - нормаль будет направлена "внутрь", к заряду. Если поверхность пересекает трубку как-то более сложно, то разобьём её на более мелкие трубки - всё равно для интегрирования нам надо будет взять предел, при котором

Далее, поставим каждой трубке

такой элемент поверхности

который высекается из поверхности интегрирования этой трубкой. Если трубка достаточно тонкая, то поверхность интегрирования в этом элементе приблизительно плоская, и может быть задана одним вектором нормали

а расстояние от точки заряда - расстоянием

и радиус-вектором

Тогда поперечное сечение трубки в месте пересечения будет

а площадь элемента поверхности - с учётом взаимного наклона поверхности и трубки,

где

- угол наклона поверхности, отсчитываемый так, что нормальная поверхность имеет

- тогда

Посчитать

можно по правилам векторной алгебры, поскольку мы знаем, что

Соберём всё это вместе:

И подставим в интеграл Гаусса для нашего одиночного точечного заряда:

У нас остался интеграл от телесных углов, взятый по всей поверхности, охватывающей заряд. Из-за того, что поверхность замкнута, и охватывает заряд ровно один раз, этот интеграл равен

даже если трубка пересекает поверхность несколько раз: при пересечениях "нормалью внутрь"

будет формально отрицательным, и в итоге вклад каждой трубки в сумме по всем пересечениям поверхности будет ровно равен угловому размеру

этой трубки. Итак:

II. Точечный заряд находится вне поверхности интегрирования.
Проделаем всё точно так же: разобьём всё пространство вокруг заряда на тонкие конические трубки. Но теперь поверхность интегрирования будет пересекать каждую трубку чётное число раз: на

пересечений "нормалью наружу" приходится столько же,

пересечений "нормалью внутрь". Далее все выкладки аналогичны, кроме того, что интеграл от телесных углов,
взятых с учётом направления нормали, будет равен 0. По сути, мы видим "из точки заряда" поверхность интегрирования снаружи, и угловая площадь "её стороны, повёрнутой к нам", будет в точности равна угловой площади "стороны, повёрнутой от нас", так как поверхность замкнутая. Итак:

Суммируем по всем зарядам, и получаем искомое соотношение.