Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to implement the bulk ferromagnetic hysteresis model (mostly the Jiles-Atherton Model), see http://drum.lib.umd.edu/bitstream/1903/6043/1/PhD_99-1.pdf page 44 equation (30).

The needed equations are $$ \delta = \text{sgn}(\dot{H})\\[1ex] H_\text{eff}=H+\alpha M \\[1ex] M_\text{an}(H_\text{eff})= M_\text{sat}\cdot\mathcal{L}(H_\text{eff}/a)\\[1ex] \delta_\text{M}=\begin{cases}0 & \dot{H}<0\text{ and }M_\text{an}-M>0\\0 & \dot{H}>0\text{ and }M_\text{an}-M<0\\1&\text{else}\end{cases}$$ and the main equation is $$ \frac{\text{d}M}{\text{d}H}=\frac{\frac{k\delta c}{\mu_0}\frac{\text{d}M_\text{an}}{\text{d}H_\text{eff}}+\delta_\text{M}(M_\text{an}-M)}{\frac{k\delta}{\mu_0}-\delta_\text{M}(M_\text{an}-M)\alpha-\frac{k\delta \alpha c}{\mu_0}\frac{\text{d}M_\text{an}}{\text{d}H_\text{eff}}}$$

Short explanation: $\delta$ takes into account wheter we are at the upper hysteresis branch or the lower one, $H_\text{eff}$ is needed to shift the Langevinfunction $\mathcal{L}(x)$ to achieve a hysteresis, $M_\text{an}$ is the anhysteretic Magnetization and $\delta_\text{M}$ is necessary to prevent unphysical behaviour.

How I tried to implement it

There were two approaches I tried: The first one is based on ParametricNDSolve, and, as I understand it, the values need to go from $0$ to $H_\text{max}$ (this is no problem and works) then down to $-H_\text{max}$ and up again to $H_\text{max}$ (this is the part I'm unsure how to achieve).

My first guess was to define the function piecewise, with an implemented shift (Beware monster ugly function below!):

\[Micro]0 = 4 \[Pi]*10^-7;
BulkFerromagneticHysteresisModelComplete = ParametricNDSolveValue[{Mges'[Happlied] == 
 Piecewise[{{((k*c/\[Micro]0)*
        Msat (a/(Happlied + \[Alpha]*Mges[Happlied])^2 - 
          Csch[(Happlied + \[Alpha]*Mges[Happlied])/a]^2/a) + 
       Piecewise[{{0, 
           Msat*Piecewise[{{#/3, 
                Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied + \[Alpha]*Mges[Happlied])/
                a) - Mges[Happlied] > 0}, {1, 
           True}}]*(Msat*
             Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied + \[Alpha]*Mges[Happlied])/
             a) - Mges[Happlied]))/((k/\[Micro]0) - 
       Piecewise[{{0, 
           Msat*Piecewise[{{#/3, 
                Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied + \[Alpha]*Mges[Happlied])/
                a) - Mges[Happlied] > 0}, {1, 
           True}}]*(Msat*
             Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied + \[Alpha]*Mges[Happlied])/
             a) - Mges[Happlied])*\[Alpha] - (k*
          c*\[Alpha]/\[Micro]0)*
        Msat (a/(Happlied + \[Alpha]*Mges[Happlied])^2 - 
          Csch[(Happlied + \[Alpha]*Mges[Happlied])/a]^2/a)), 
    Happlied < 
     20000}, {((-k*c/\[Micro]0)*
        Msat (a/(-Happlied + 
            40000 + \[Alpha]*Mges[-Happlied + 40000])^2 - 
          Csch[(-Happlied + 
            40000 + \[Alpha]*Mges[-Happlied + 40000])/a]^2/a) + 
       Piecewise[{{0, 
           Msat*Piecewise[{{#/3, 
                Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((-Happlied + 
                40000 + \[Alpha]*Mges[-Happlied + 40000])/a) - 
             Mges[-Happlied + 40000] > 0}, {1, 
           True}}]*(Msat*
             Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((-Happlied + 
               40000 + \[Alpha]*Mges[-Happlied + 40000])/a) - 
          Mges[-Happlied + 40000]))/((-k/\[Micro]0) - 
       Piecewise[{{0, 
           Msat*Piecewise[{{#/3, 
                Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((-Happlied + 
                40000 + \[Alpha]*Mges[-Happlied + 40000])/a) - 
             Mges[-Happlied + 40000] > 0}, {1, 
           True}}]*(Msat*
             Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((-Happlied + 
               40000 + \[Alpha]*Mges[-Happlied + 40000])/a) - 
          Mges[-Happlied + 40000])*\[Alpha] - (-k*
          c*\[Alpha]/\[Micro]0)*

        Msat (a/(-Happlied + 
            40000 + \[Alpha]*Mges[-Happlied + 40000])^2 - 
          Csch[(-Happlied + 
            40000 + \[Alpha]*Mges[-Happlied + 40000])/a]^2/a)), 
    20000 < Happlied < 60000},{((k*c/\[Micro]0)*
        Msat (a/(Happlied + 
            40000 + \[Alpha]*Mges[Happlied - 80000])^2 - 
          Csch[(Happlied + 
            40000 + \[Alpha]*Mges[Happlied - 80000])/a]^2/a) + 
       Piecewise[{{0, 
           Msat*Piecewise[{{#/3, 
                Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied -80000 + \[Alpha]*Mges[Happlied -80000])/a) - 
             Mges[Happlied -80000] > 0}, {1, 
           True}}]*(Msat*
             Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied -80000 + \[Alpha]*Mges[Happlied -80000])/a) - 
          Mges[Happlied -80000]))/((k/\[Micro]0) - 
       Piecewise[{{0, 
           Msat*Piecewise[{{#/3, 
                Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied - 
                80000 + \[Alpha]*Mges[Happlied - 80000])/a) - 
             Mges[Happlied - 80000] > 0}, {1, 
           True}}]*(Msat*
             Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), 
                True}}] &@((Happlied - 
               80000 + \[Alpha]*Mges[Happlied -80000])/a) - 
          Mges[Happlied -80000])*\[Alpha] - (k*
          c*\[Alpha]/\[Micro]0)*

        Msat (a/(Happlied - 
            80000 + \[Alpha]*Mges[Happlied - 80000])^2 - 
          Csch[(Happlied - 
            80000 + \[Alpha]*Mges[Happlied -80000])/a]^2/a)), 
    60000 < Happlied < 100000}}], 
Mges[0] == 0.0001}, {Mges}, {Happlied, 0, 100000}, {Msat, a, 
k, \[Alpha], c}];

Unfortunately this doesn't work (I suspect because of Mges[Happlied-80000] respectively Mges[-Happlied+40000]? If yes, how can I workauround this?).

EDIT: To be more precise, compiling the above shown works, proceeding with

modelcomplete[Msat_, a_, k_, \[Alpha]_, c_][Happlied_] := Through[BulkFerromagneticHysteresisModelComplete[Msat, a, 
  k, \[Alpha], c][Happlied], List] /; And @@ NumericQ /@ {Msat, a, k, \[Alpha], c, Happlied};

leads after several seconds (it looks as if it is already done for some time with compiling) to the error:

Mathematica Kernel for Windows has stopped working
Problem Event Name: APPCRASH
Application Name: MathKernel.exe
...
Fault Module Name: mathdll.dll
...

Another possible approach is described in Computer-aided Analysis of Electric Machines: A Mathematica Approach from Vlado Ostovic, but the book is from the early 90s, so the code presented there doesn't work - the main idea is to use a periodic function to achieve the hyteresis effect. It seems I got the virgin curve working there as well:

Ms = 2*10^6;
a = 1100;
kon = 0.0004;
alpha = 1.6*10^3;
mu0 = N[4*Pi*10^(-7)];
c = 0.2;
Hmax = 8000;
omega = N[2*Pi];

eq = {Mges'[t] == (Ms (Coth[(Hmax*Cos[2 Pi t] + alpha Mges[t])/a] - 
      a/(Hmax*Cos[2 Pi t] + alpha Mges[t])) - 
   Mges[t])/((1 + c) (kon/mu0 - 
     alpha (Ms (Coth[(t + alpha Mges[t])/a] - 
           a/(Hmax*Cos[2 Pi t] + alpha Mges[t])) - Mges[t]))) + 
c/(1 + c) (Ms (a/(Hmax*Cos[2 Pi t] + alpha Mges[t])^2 - 
     Csch[(Hmax*Cos[2 Pi t] + alpha Mges[t])/a]^2/a)), Mges[0] == 0}

tests = NDSolve[eq, {Mges}, {t, 0, 1}]

Plot[Evaluate[Mges[t] /. tests], {t, 0, 1}, PlotRange -> All]

But instead of the Mges[t] it would be necessary to use Mges[Hmax*Cos[2 Pi t]], sadly Mathematica crashes if I try to use this expression.

EDIT: To be more precise here: If I use Mges[Hmax*Cos[2 Pi t]] throughout eq I get the following error:

NDSolve::derivs: No derivatives of dependent variables were found in the equations. NDSolve is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations. >>

I suspect this is due to Mathematica not recognizing Mges[Hmax*Cos[2 Pi t]] as $M(x)$ but as $M(f(x))$.


Does anyone have any idea how to overcome these obstacles? Thanks a lot in advance, I really appreciate any help!

share|improve this question
    
John, could you please clarify what you mean by “doesn’t work” for the first method? And for the second method do you mean the kernel crashes? –  Verbeia Apr 26 at 16:56
    
@Verbeia I tried to clarify the admittedly unclear formulation, see the two edits. Thanks for pointing out! –  John Apr 26 at 17:19

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.