Модель гистерезиса Прейзаха

TL;DR

В конце этого длинного текста описано создание физического компонента на языке физического моделирования Engee. Можно сразу попробовать помоделировать цепи с этим компонентом. Или почитать зачем и что именно сделано, начиная с мотивации и быстрой теории.

Мотивация

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

Это связано со сложностью самого явления, его описания и, при этом, незначительными улучшениями результатов. Периодически мелькают статьи, где различия в форме вторичного тока трансформатора тока при учете гистерезиса и без учета невелики.

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

1. Остаточная намагниченность

   С однозначной зависимостью между током и потоком невозможно получить результат, когда в начале процесса при нулевом токе вы имеете нулевой поток, а в конце процесса при нулевом токе уже остается какой-то поток. В этом суть однозначной зависимости - при определенном значении тока жестко определенный поток. Остаточную намагниченность можно получить только с учетом гистерезиса.

   Это важно потому, что современные ГОСТы по расчету времени до насыщения ТТ предлагают учитывать остаточную намагниченность в 86% (для замкнутого магнитопровода). Это довольно высокое значение, которое сильно влияет на результат расчета. Такой уровень остаточной намагниченности периодически критикуется сообществом. Верифицированная модель ТТ с учетом гистерезиса позволила бы в рамках численных экспериментов набрать статистику по остаточной намагниченности (фаза КЗ хоть и случайная величина, но распределена неравномерно).

2. Размагничивание ТТ токами нагрузки

   После возникновения КЗ и его отключения, ТТ нужно размагничивать, чтобы он вернулся в класс точности. Тем не менее в статьях [1] указывается, что размагничивание происходит и токами нагрузки, причем именно при длительном их воздействии. Объяснение этого эффекта невозможно даже при учете гистерезиса стандартными моделями. Нужны более совершенные модели.

Базовые понятия

Вспомним базовые вещи про петлю гистерезиса. Внешне она может выглядеть примерно так:

Это известные всем параметры.

Так же часто полезна величина . Это значение напряженности поля выше которого (по модулю) гистерезисные свойства уже не проявляются, индукция растет медленно и можно считать, что зависимость становится однозначной. Также при возврате напряженности поля из области с большими перемагничивание будет идти по нисходящей (восходящей) ветви предельной петли гистерезиса.

Кривые возврата

Если во время движения по предельной нисходящей ветви (при уменьшении напряженности поля) поменять направление намагничивания (начать увеличивать напряженность поля) то перемагничивание будет осуществляться по кривой возврата первого порядка (FORC - first order reversal curve). Так называются кривые которые начинаются на предельной петле. Можно покрыть петлю гистерезиса двумя семействами таких кривых, начинающихся на восходящей и нисходящей ветви предельной петли. Кривые, которые начинаются на FORC называются кривыми возврата второго порядка и так далее. Модели гистерезиса можно классифицировать по "умению" моделировать кривые возврата и по необходимости для идентификации модели иметь набор таких экспериментальных кривых.

Правила Маделунга (von E. Madelung)

Обобщение экспериментальных данных о том, как перемагничивается материал при различных входных воздействиях, принципы построения частных циклов перемагничивания были изложены в статье [2] в виде четырех правил. Трактовались они не всегда однозначно. Я их изложу согласно книге [3] и статье [4].

Правило 1

Ход частной кривой B(H) определяется в общем случае всей предысторией процесса перемагничивания.

Это некоторое базовое утверждение, по сути, еще одно определение явления гистерезиса.

Правило 2

Если некоторую точку 3 кривой возврата 2-3-1, начинающейся в точке поворота 2, сделать новой точкой поворота, то кривая 3-4-2 возвращается в точку 2 предыдущего поворота.

Правило 3

Если точка 4 кривой 3-4-2 в свою очередь становится точкой поворота, то кривая 4-3, пройдя через точку 3, совпадёт при своем продолжении с отрезком 3-1 предшествующей кривой 2-3-1, как если бы замкнувшегося цикла 3-4-3 вовсе не существовало.

Правила 2 и 3 качественно определяют, как должны строиться кривые перемагничивания. Правило 4, как наиболее приближенное я здесь указывать не буду. Его выполнение / не выполнение может меняться от выбранной модели гистерезиса.

Свойства известных моделей

Модель Джайлса-Атертона (Jiles–Atherton model)

Когда говорят про гистерезис часто подразумевают модель Джайлса-Атертона. Она реализована в программных средах и, как будто, наиболее популярна.
С другой стороны у нее есть ряд минусов и она обрастает модификациями. Я не являюсь специалистом по этой модели и не исключаю, что нужно уметь ее настраивать. Просто ограничусь некоторыми примерами из статей.

В статье [5] указаны возможные отрицательные наклоны кривой намагничивания, что нефизично.

В статье [6] указано на возможное самопересечение частного цикла и предложена соответствующая модификация. Такое самопересечение нарушает 2-ое правило Маделунга.

В ряде источников [7] указана неудовлетворительная форма частных циклов перемагничивания.

Также не самым простым является подбор параметров модели под экспериментальные данные.

Модель Джона Чана (John H. Chan)

В данной модели и в ее модификациях [8] не выполняется правило 2 о замыкании частного цикла. При развороте в точке D и последующем развороте в точке F кривая должна вернуться в точку D, а не прийти в новую точку G.

Модель гистерезиса Прейзаха лишена всех описанных минусов (почти) и “из коробки” соответствует правилам Маделунга. К ней и переходим.

Модель гистерезиса Прейзаха

Здесь я дам краткое описание, чтобы потом было проще перейти к описанию кода физического компонента. Отличное изложение теории этой модели дано в книге [9].

Базовым элементом модели является гистерон. Это элемент который включается (+1), когда напряженность поля превышает порог включения и отключается (-1), когда поле снижается ниже порога отключения . Каждый гистерон имеет собственный вес. Гистерон - элементарное магнитное реле.

Магнитный материал представялет собой множество таких гистеронов со всевозможными порогами включения/отключения и специально подобранными весами. Суммарное действие всех гистеронов определяет как величину поля (индукция), так и сложное нелинейное поведение (память) при изменении напряженности поля.

Треугольник Прейзаха

Рассмотрим геометрическую интерпретацию модели Прейзаха. Это основная абстракция, которая помогает понять свойства модели и численно её реализовать.

Будем располагать гистероны с порогами и на плоскости. По оси X будем располагать пороги отключения , по оси Y пороги включения . Порог включения всегда выше порога отключения то есть > . Это неравенство определяет полуполскость в координатах , . Так же, как мы ввели ранее, есть величина выше которой зависимость становится однозначной. То есть, если напряженность поля выше то все гистероны включены, если ниже минус , то все гистероны выключены. Этими значениями ограничивается треугольник.

Каждая точка внутри треугольника это какой то гистерон. Введем функцию на этой плоскости . Значение функции в конкретной точке показывает вес гистерона в этой точке.
Далее я буду часто рисовать этот треугольник. Включенные гистероны я буду обозначать жёлтым цветом, выключенные синим.

Как работать с треугольником Прейзаха

Увеличим напряженность поля, например, до . Магнитный материал полностью насытится и значение индукции будет равно . Все гистероны включатся, так как превышает пороги включения всех гистеронов. Значение должно быть сформировано как сумма весов всех гистеронов (в непрерывном случае значению должен быть равен двойной интеграл от по области, ограниченной треугольником).
Мы получим такую картинку.

Теперь уменьшим напряженность поля до минус . Магнитный материал также насытится, но с другим знаком. Все гистероны выключатся. Их суммарный вклад даст минус .
Картинка будет такой.

Теперь увеличим поле до . Включатся все гистероны для которых > . Геометрически это означает, что включатся все, которые лежат ниже горизонтальной прямой с ординатой . В координатах B(H) мы пройдем от отрицательного насыщения до по восходящей предельной кривой намагничивания.

Уменьшим поле до минус . При уменьшении поля из включенных гистеронов начинают отключаться только те, для которых текущая напряженность поля становится меньше их порога возврата . Геометрически это означает, что выключаются гистероны правее вертикальной прямой с абсциссой минус . Важно, что отключается только часть гистеронов (в силу индивидуального гистерезиса) из тех, что были включены. А значит при уменьшении напряженности поля вклад в уменьшение магнитной индукции будет меньше, чем при нарастании напряженности поля и в координатах мы увидим движение по кривой возврата (FORC), которая отличается от восходящей кривой.

Опять поменяем направление изменения напряженности поля. Надеюсь картинки уже сразу понятны.

Тут самым важным является то, что по мере движения горизонтальной прямой вверх (при нарастании поля) рано или поздно мы совпадем с предыдущей "горизонтальной" ступенькой. В этот момент материал попадает в состяние намагниченности в котором уже был. А значит кривая в координатах пришла ровно в предыдущую точку поворота (Правило 2 выполняется). При дальнейшем увеличении напряженности поля двигать уже надо увеличенную "ступеньку". Это означает, что точка будет двигаться ровно по той же траектории что и до формирования этого частного цикла (Правило 3 выполняется).

Таким образом магнитное состояние образца полностью определяется лестничного вида ломаной линией, разделяющей треугольник на включенные и отключенные гистероны. А "слияние" вертикальных или горизонтальных частей ступеней при изменении напряженности поля позволяет сразу выполнить правила Маделунга. Этим данная модель и отличается от остальных.

Минусы классической модели Прейзаха

Главным минусом данной модели является то, что она порождает конгруэтные частные циклы перемагничивания которые ограничены одними и теми же значениями напряженности поля. То есть любой замкнутый цикл можно сдвигать вверх или вниз и он будет оставаться корректным. Это приводит к тому, что замечание по неудовлетворительной форме частных циклов относится и к данной модели. Решение этого вопроса может быть получено в более сложной нелинейной модели Прейзаха.

Идентификация функции Прейзаха

Постоянное распределение

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

Получается вот такой параболических гистерезис, где любая кривая перемагничивания от одной точки поворота до другой - парабола с вершиной в точке разворота.

Кусочно постоянное распределение

Теперь рассмотрим распределение гистеронов, где на треугольнике справа и снизу гистероны имеют нулевой вес (или можно считать что их там нет вообще). А гистероны в оставшемся квадрате имеют постоянный вес.

При таком распределении движение к нулевой напряженности поля (по оси снизу вверх и по оси справа налево) никак не будет менять магнитную индукцию. Следовательно в осях это горизонтальные линии. При движении через квадрат (или часть квадрта) в силу одинаковой высоты площадь (и объем) будет всегда меняться линейно. Получается своеобразный угловатый гистерезис, в котором тем не менее продолжают выполняться правила Маделунга.

Идентификация общего вида

Функция распределения гистеронов это ядро модели Прейзаха. Поэтому ставится вопрос о выборе (расчете) функции распределения гистеронов (или в дискретном виде значений на квадратной сетке), которая бы описывала реальные данные. Возможность однозначной идентификации зависит от количества исходных данных. Если даны только предельные петли, то однозначно задача не решается и возможен выбор из множества подходящих функций. В физическом компоненте как раз реализован такой вариант. Идентификация матрицы идёт только по параметрам предельной петли. Про функцию распределения стоит добавить, что она должна быть всегда неотрицательной и для порождения симметричных относительно начала координат предельных циклов сама являться симметричной .

Один из вариантов реальной идентификации показан ниже.

Физический компонент

Реализуем вышеописанную модель как физический компонент, с помощью языка физического моделирования Engee.
Нам необходимо держать в памяти всю ступенчатую линию, которая разделяет треугольник Прейзаза, так как именно она совместно с функцией распределения определяет текущую величину потока (индукции). Также необходимо реализовать логику добавления и удаления ступенек при увеличении/уменьшении напряженности магнитного поля, когда мы выходим за рамки текущей ступеньки.

Лого

Но начнем мы с самого главного в физическом компоненте, конечно с картинки. Подберем нужный цвет (чтобы попасть в магнитный домен Engee) и разработаем лого. Исходя из большого количества треугольников, которые я нарисовал выше, считаю, что логичным обозначением гистерезиса Прейзаха будет всё таже полуплоскость.

Код

После объявления физического компонента как такового через @engeemodel объявляем ненаправленные порты. Обозначим их как северный и южный:

@nodes begin
    p = EngeePhysicalFoundation.Magnetic.Port,
    [label="N",
    side="left",
    offset=-7]

    n = EngeePhysicalFoundation.Magnetic.Port,
    [label="S",
    side="right",
    offset=-7]
end

Главными переменными будут поток и магнитное напряжение (магнитодвижущая сила):

@variables begin
    mmf = 0.0, [unit = "A", description = "Магнитодвижущая сила"]
    flux = 0.0, [unit = "Wb", description = "Поток", connect = Flow]
end

Введем промежуточную переменную (производную от МДС), чтобы определять максимумы (минимумы) на кривой МДС и добавлять ступеньки в ломаную ступенчатую линию, разделяющую треугольник Прейзаха:

@intermediates begin
    der_mmf = der(mmf)
end

Определим первые параметры. Для задания гистерезса требуется задать линию (значения по x и y) для восходящей предельной ветви и наклон в однозначной части. Матрицу распределения весов компонент вычислит сам.

@parameters [group = "Данные по восходящей предельной петле намагничивания"] begin
    mmf_data[:] = [-1.0, -0.925246908773959, -0.825766786751853, -0.70872836141771, -0.581300360255558, -0.450651510749425, -0.323950540383338, -0.208366176641324, -0.111067147007411, -0.0392221789656274, 0.0, 0.0186377772808986, 0.0338781632394863, 0.0462753648482215, 0.0563835890795626, 0.0647570429059679, 0.0719499332998958, 0.0785164672338045, 0.0850108516801526, 0.0919872936113983, 0.1, 0.109845234455496, 0.125060359224361, 0.150311229234456, 0.190263699413644, 0.249583624689787, 0.332936859990747, 0.444989260244386, 0.590406680378567, 0.773854975321151, 1.0], [description="Магнитное напряжение, А"]
    flux_data[:] = [-1.0, -0.999474340152244, -0.997242098319089, -0.992320381065707, -0.983726294957268, -0.970476946558943, -0.951589442435903, -0.926080889153317, -0.892968393276358, -0.851269061370195, -0.8, -0.740896690912631, -0.676509837062578, -0.607079376826463, -0.532845248580905, -0.454047390702525, -0.370925741567943, -0.283720239553779, -0.192670823036654, -0.0980174303931872, 0.0, 0.106730766333609, 0.224540461935931, 0.348713591136731, 0.474534658265773, 0.597288167652825, 0.712258623627651, 0.814730530520016, 0.899988392659685, 0.963316714376425, 1.0], [description="Поток, Вб"]
    Ls = 0.001, [description="Наклон в однозначной части хар-ки"]
end

Определим количество гистеронов на катете треугольника Прейзаха (я не рекомендую ставить больше 100, в полном треугольнике будет ~5 000 штук). Также создадим выпадающий список, для задания начального состяния гистеронов:

@structural_parameters [group = "Данные по модели Прейзаха"] begin
    n_hysterons::Int64 = 50, [description="Количество гистеронов на катете треугольника Прейзаха"]
    choosenPoint::startPoint = startPoint.neg, [description="Выбор нач. состояния"]
end

Для того чтобы выпадающий список работал, до начала кода компонента нужно вставить:

@descriptive_enumeration startPoint begin
    neg = 1, "Отрицательное насыщение"
    pos = 2, "Положительное насыщение"
    zero = 3, "Размагниченное состояние"
end

Следующий кусок кода (условие) сработает один раз (при инициализации) и, в зависимости от выбранного начального состояния, немного по-разному определит дополнительные параметры и событийные переменные.

if choosenPoint == startPoint.neg
    @parameters [gui = None, access = Private] begin
        extrems[:] = [mmf_sat, -mmf_sat]
    end

    @intermediates [event=true] begin
        increase = true
        last_min = -mmf_sat
        last_max = mmf_sat
    end

elseif choosenPoint == startPoint.pos
    @parameters [gui = None, access = Private] begin
        extrems[:] = [mmf_sat, -mmf_sat, mmf_sat]
    end

    @intermediates [event=true] begin
        increase = false
        last_min = -mmf_sat
        last_max = mmf_sat
    end

else
    @parameters [gui = None, access = Private] begin
        extrems[:] = demagnetization(mmf_sat)
    end

    @intermediates [event=true] begin
        increase = true
        last_min = extrems[end]
        last_max = extrems[end-1]
    end
end

Параметр extrems это динамический массив, в котором будут запоминаться вершины ломаной ступенчатой линии. В него будем добавлять значения при формировании новой точки поворота кривой намагничивания и удалять значения при "затирании" старых ступенек.
Булевская переменная increase будет использоваться для реализации логики управления массивом.

Дополнительные скрытые параметры. Остановиться здесь можно на матрице discreteTriangle, в которой записаны значения весов гистеронов (функция распределения). Матрица определяется как результат выполнения функции identify_matrix(). Сама функция определена перед кодом компонента:

@parameters [gui = None, access = Private] begin
    mmf_sat = mmf_data[end]
    flux_sat = flux_data[end]
    discreteTriangle[:, :] = identify_matrix(n_hysterons, mmf_data, flux_data)
    last_extrems[:] = [mmf_sat, -mmf_sat]
end

Стандартный блок - аналог первого закона Кирхгофа для магнитных цепей. Здесь отмечу, что в хелпе показаны только электрические компоненты и то, что через точку обращаться к свойству порта нужно именно через имя Phi сразу непонятно.

@branches begin
    flux:(p.Phi, n.Phi)
end

Блок с уравнениями. Условия в этом блоке проверяются на каждом шаге расчета физической модели.

@equations begin
    mmf ~ p.mmf - n.mmf

    if mmf >= mmf_sat
        flux ~ flux_sat + Ls * (mmf - mmf_sat)
    elseif mmf <= -mmf_sat
        flux ~ -flux_sat + Ls * (mmf + mmf_sat)
    else
        flux ~ Symbolics.term(get_flux, extrems, mmf, mmf_sat, flux_sat, discreteTriangle, increase, everett_fun, type=Real)
    end
end

Весь гистерезис спрятан в функцию get_flux(), которая как и другие фукнции определена до компонента и вызывается через специальную обертку Symbolics.term(). Она нужна, чтобы символьный движок не пытался залезть внутрь и разложить ее на более элементарные части.

Событийные переменные

События позволяют дискретно (в момент события) менять поведение компонента через уравнения и параметры. Через блоки @when можно фиксировать события. Пойдем по порядку:

1. Если МДС превысила МДС насыщения (переход на однозначную часть зависимости), то сбрасываем значения массива extrems функцией pos_reset!(). Принципиальным тут является тот факт, что мы меняем вектор extrems, который является параметром модели, а не переменной. Переменная же dummy нужна просто для того, чтобы функция могла куда-то вернуть скалярное значение (выглядит как костыль, но так работает);
2. Второе условие - тоже самое, но для отрицательного насыщения;
3. Если значение МДС стало выше последней ступеньки, то ступеньки надо "слить" в одну. Этим занимается функция del_pos_extremum!(). Её главная задача также преобразовывать extrems;
4. Четвертое условие - тоже самое при уменьшении МДС;
5. Если производная МДС поменяла знак - значит в процессе намагничивания появилась точка разворота, то нужно добавлять ступеньку в extrems функцией save_extremmum!();
6. Тоже самое для изменения знака производной в другую сторону.

@when edge(mmf >= mmf_sat) | edge(t > 0.0) & ((mmf >= mmf_sat) & (increase == true) | (mmf > last_extrems[1]) & (increase == false)) begin
    increase ~ false
    dummy ~ Symbolics.term(pos_reset!, extrems, mmf_sat, last_extrems, type=Real)

@elsewhen edge(mmf <= -mmf_sat) | edge(t > 0.0) & ((mmf <= -mmf_sat) & (increase == false) | (mmf < last_extrems[2]) & (increase == true))
    increase ~ true
    dummy ~ Symbolics.term(neg_reset!, extrems, mmf_sat, last_extrems, type=Real)

@elsewhen edge(mmf > (last_extrems[1]+0.001)) & (increase == true)
    # Обрезаем массив extrems
    dummy ~ Symbolics.term(del_pos_extremum!, extrems, mmf, last_extrems, type=Real)

@elsewhen edge(mmf < (last_extrems[2]-0.001)) & (increase == false)
    # Обрезаем массив extrems
    dummy ~ Symbolics.term(del_neg_extremum!, extrems, mmf, last_extrems, type=Real)

@elsewhen edge(der_mmf > 0.0) & (increase == false) | edge(t > 0.0) & (der_mmf > 0.0) & (increase == false)
    # Нашли локальный минимум
    increase ~ true
    dummy ~ Symbolics.term(save_extremum!, extrems, mmf, last_extrems, type=Real)

@elsewhen edge(der_mmf < 0.0) & (increase == true) | edge(t > 0.0) & (der_mmf < 0.0) & (increase == true)
    # Нашли локальный максимум
    increase ~ false
    dummy ~ Symbolics.term(save_extremum!, extrems, mmf, last_extrems, type=Real)
end

На первом шаге, решатель ничего не знает про производную МДС, поэтому необходимо добавить:

@defaults begin
    der(mmf) => 1e-6
end

Не все вышеприведенные конструкции описаны в хелпе. Без помощи в чате техподдержки ничего бы не получилось. Спасибо за неё.

Как помочь вычислителю?

Численное решение системы уравнений, описывающей физическую систему подразумевает на каждом шаге решение нелинейной системы уравнений. Базой для решения является метод Ньютона, следовательно нужна матрица Якоби. В контексте данного компонента нужна производная потока по магнитному напряжению (дифференциальная индуктивность, дифференциальная магнитная проницаемость). Без специального учета этого фактора алгоритм на этапе инициализации будет пытаться определить производную самостоятельно, залезая в функцию get_flux(), что у меня всегда приводило к ошибкам. Поправить ситуацию можно введением до кода компонента следующих строк:

Symbolics.@register_symbolic der_flux_mmf(mmf, extrems, M, increase, mmf_sat)
Symbolics.derivative(::typeof(get_flux), args::NTuple{7,Any}, ::Val{2}) = Symbolics.term(der_flux_mmf, args[2], args[1], args[5], args[6], args[3])

В первой строке мы регистрируем функцию der_flux_mmf(), которая по треугольнику Прейзаха вычисляет производную потока (практически аналитическое задание производной). Регистрируем её для того, чтобы она рассматривалась, как базовая функция и алгоритм не пытался разложить ее на более элементарные куски. Во второй строке мы определяем производную функции get_flux() по второму аргументу (это как раз mmf). То есть теперь, когда алгоритму понадобится производная, он просто вызовет соответствуюущую функцию и не будет пытаться самостоятельно ее определить.

Применение магнитного компонента

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

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

Электрический компонент

Признаемся что популярность магнитного домена оставляет желать лучшего. Из 144-х 1-D моделей в сообществе Engee магнитный домен использован в парочке штук. Поэтому на базе реализованного гистерезиса разработаем уже электрическую катушку с гистерезисом. Посмотрим только на отличия в коде от магнитного компонента.

Лого

Во-первых, конечно, разработаем другой логотип:

Код

Определяем другие порты:

@nodes begin
    p = EngeePhysicalFoundation.Electrical.Pin,
    [label="+",
    side="left",
    offset=-7]

    n = EngeePhysicalFoundation.Electrical.Pin,
    [label="-",
    side="right",
    offset=-7]
end

Другие переменные (электрические, а не магнитные)

@variables begin
    v = 0.0, [unit = "V", description = "Напряжение"]
    i = 0.0, [unit = "A", description = "Ток", connect = Flow]
end

Поток flux теперь переезжает в промежуточные переменные. Также определяем mmf с учётом количества витков катушки w.

@intermediates begin
    der_i = der(i)
    mmf = i*w
    flux = ifelse(mmf >= mmf_sat, flux_sat + Ls * (mmf - mmf_sat), ifelse(mmf <= -mmf_sat, -flux_sat + Ls * (mmf + mmf_sat), Symbolics.term(get_flux, extrems, mmf, mmf_sat, flux_sat, discreteTriangle, increase, everett_calc2, type=Real) ))
end

Меняем переменные при описании ветви:

@branches begin
    i:(p.i, n.i)
end

И самые главные изменения в секции уравнений:

@equations begin
    v ~ p.v - n.v
    v ~ ifelse(mmf >= mmf_sat, w * Ls * der_i, ifelse(mmf <= -mmf_sat, w * Ls * der_i, w * der_i * Symbolics.term(der_flux_mmf2, mmf, extrems, discreteTriangle, increase, mmf_sat, type=Real)))
end

Уравнение для напряжения имеет вид . Но если его оставить в таком виде, то вычислитель будет опять делать попытки пролезть в функцию get_flux(). Если же мы применим правило дифференцирования сложной функции, то получим нелинейное дифференциальное уравнение уже на ток, тем более производную потока по току (МДС) мы уже вычисляли ранее.

В таком виде компонент уже рассчитывается.

Применение электрического компонента

Применим компонент в такой схеме:

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

Заключение

С помощью языка физического моделирования разработано два компонента, реализующих классическую модель гистерезиса Прейзаха в магнитном и электрическом доменах. Компоненты работают, но сломать расчет всё еще легко. Часто ошибки лечатся рекомендацией из хелпа, но не всегда. Не хватило режима отладки, когда в моменте физического моделирования можно легко посмотреть на любую переменную компонента, в том числе внутри отдельно объявленных функций.

Литература

[1] Определение остаточных магнитных индкуций в тороидальных сердечниках трансформаторов тока класса Р для релейной защиты. Дегтяров А.А., Кужеков С.Л., Дони Н.А., Шурупов А.А.
[2] Uher Maynetdsderung durch schmellverlaufende StrSme und dde Wdrkunyswedse des Rutherford-Narc owdschen Magnetdetektors. Von E. Madelung
[3] Переходные режимы работы ТТ. Сирота И.М.
[4] Моделирование магнитного гистерезиса на основе обобщенных правил Маделунга. Часть 1. Зирка С.Е., Мороз Ю.И.
[5] A differential equation approach to minor loops in the Jiles-Atherton Hysteresis Model. Carpenter, 1991
[6] Minor loops modelling with a modified Jiles–Atherton model and comparison with the Preisach model. Benabou, 2008
[7] Micro-Cap. Версии 9 и 10. Амелин С.А., Амелина М.А.
[8] Гистерезисная модель нелинейной индуктивности симулятора LTspice. Володин Валентин.
[9] Mathematical models of hysteresis and their applications. I. Mayergoyz.