Archived OpenModelica forums. Posting is disabled.

Alternative forums include GitHub discussions or StackOverflow (make sure to read the Stack Overflow rules; you need to have well-formed questions)


Forgot password? | Forgot username? | Register

Jiles Atherton hysteretic model equation problem

Jiles Atherton hysteretic model equation problem

Hello,
I am new in Modelica and Ihave problem with Jiles Atherton hysteretic model. Here is my code:

model Hysteretic_curve
  parameter Real N = 1;
  parameter Real l = 0.816;
  parameter Real Ms = 1500000.0;
  parameter Real a = 350;
  parameter Real alpha = 0.0007;
  parameter Real c = 0.001;
  parameter Real k = 265;
  constant Real pi = 3.14;
  discrete Real i(start = 0),H(start = 0),delta(start = 0),Man(start = 0),M(start = 0),B(start=0);
equation
  i = sqrt(2) * 2000 * sin(314 * time);
  H = N * i / l;
  Man = Ms * (1 / tanh((H + alpha * M) / a) - a / (H + alpha * M));
  if der(H) > 0 then
    delta = 1;
  else
    delta = -1;
  end if;
  der(M) = ((1 - c) * der(H) * (Man - M)) / ((1 - c) * k - alpha * (Man - M)) + c * der(Man);
  B=1.2566e-6*(H+M);
end Hysteretic_curve;

It´s writen out 4 errors :

Translation    17:41:43        0:0-0:0    Internal error Transformation Module failed!
Translation    17:41:43        0:0-0:0    Internal error sorting equations(strongComponents failed)
Translation    17:41:43        0:0-0:0    Internal error BackendDAETransform.analyseStrongComponentBlock failed
Translation    17:41:43        0:0-0:0    Internal error Sorry - Support for Discrete Equation Systems is not yed implemented

It would be great if somebody could help. Thanks.

Re: Jiles Atherton hysteretic model equation problem

I have a practical but probably not the best solution to your problem :

Code :

model Hysteretic_curve
  parameter Real N = 1;
  parameter Real l = 0.816;
  parameter Real Ms = 1500000.0;
  parameter Real a = 350;
  parameter Real alpha = 0.0007;
  parameter Real c = 0.001;
  parameter Real k = 265;
  constant Real pi = 3.14;
  Real deridt, i, delta, B;
  Real H(start = 0),Man(start = 0),M(start = 0);
equation
  i = sqrt(2) * 2000 * sin(314 * time);
  deridt = sqrt(2) * 2000 * 314 * cos(314 * time);
  der(H) = N * deridt / l;
  //---------------------
  Man = Ms * (1 / tanh((H + alpha * M) / a) - a / (H + alpha * M));
  // Compute der(Man) with der(H) and der(M)
  // You may use Maxima to compute litteral expression
  //---------------------
  if der(H) > 0 then
    delta = 1;
  else
    delta = -1;
  end if;
  der(M) = ((1 - c) * der(H) * (Man - M)) / ((1 - c) * k - alpha * (Man - M)) + c * der(Man);
  B=1.2566e-6*(H+M);
end Hysteretic_curve;

Re: Jiles Atherton hysteretic model equation problem

Thank you! Now there are not errors, but the current is no more sunisoidal and all results are  wrong. Have you an idea, what´s wrong now?

Re: Jiles Atherton hysteretic model equation problem

Can send me
   1/ the modified code
   2/ the version of OM you use
   3/ the selected solver

I tried to calculate der(Man) and put it in the code but the my version
of OM crashed after few time steps.

Re: Jiles Atherton hysteretic model equation problem

For information, I simulate the model with the following code :

function sech
    input Real x;
    output Real y;
algorithm
    y := 2.0 / (exp(x) + exp(-x));
end sech;

model Hysteretic_curve
    parameter Real N = 1;
    parameter Real l = 0.816;
    parameter Real Ms = 1500000.0;
    parameter Real a = 350.0;
    parameter Real alpha = 0.0007;
    parameter Real c = 0.001;
    parameter Real k = 265.0;
    Real i , der_i_dt;
    Real H , der_H_dt;
    Real H_alpha_M, der_H_alpha_M_dt ;
    Real delta ;
    Real Man, der_Man_dH_alpha_M, der_Man_dt ;
    Real B ;
    Real M(start = 1.0e-6);
    Real coeff_debug1;
equation
    i        =         sqrt(2.0) * 2000.0 * sin(314.0 * time);
    der_i_dt = 314.0 * sqrt(2.0) * 2000.0 * cos(314.0 * time);
    H        = N * i / l ;
    der_H_dt = (N/l) * der_i_dt ;

    H_alpha_M        = ( H        + alpha * M      ) / a;
    der_H_alpha_M_dt = ( der_H_dt + alpha * der(M) ) / a;

    Man                = Ms * (   1.0 / tanh(H_alpha_M)                    - 1.0 /  H_alpha_M     );
    der_Man_dH_alpha_M = Ms * ( -(sech(H_alpha_M))^2 / (tanh(H_alpha_M))^2 + 1.0 / (H_alpha_M^2) ) ;
    der_Man_dt         = der_Man_dH_alpha_M * der_H_alpha_M_dt ;

    if der_H_dt > 0 then
        delta = 1;
    else
        delta = -1;
    end if;

    coeff_debug1 = (1.0 - c) * k - alpha * (Man - M) ;
    der(M) = ((1.0 - c) * der_H_dt * (Man - M)) / ((1.0 - c) * k - alpha * (Man - M)) + c * der_Man_dt;
    B      = 1.2566e-6  * ( H + M );
end Hysteretic_curve;


I think the the problem is due to coeff_debug1 wich tends to 0 during the first time steps.
The derivative of M tends to infinity and the solver is not able to compute a solution.
I used 0.000001 s time step and the crash occurs at t= 0.009021104s

The ODE of must be different if (1.0 - c) * k - alpha * (Man - M) tends to 0

Re: Jiles Atherton hysteretic model equation problem

I use OM 1.8.1 and here is my code:

model Hysteretic_curve
    parameter Real N = 1;
    parameter Real l = 0.816;
    parameter Real Ms = -1500000;
    parameter Real a = 350;
    parameter Real alpha = 0.0007;
    parameter Real c = 0.001;
    parameter Real k = 265;
    Real deridt,i,delta,B,derH,derMan;
    Real H(start = 0),Man(start = 0),M(start = 1e-005);
equation
    i = 2000 * sin(314 * time);
    deridt = 2000 * 314 * cos(314 * time);
    derH = (N * deridt) / l;
    der(H) = derH;
    derMan = Ms * ((derH + alpha * der(M)) / (a * (sinh((H + alpha * M) / a)) ^ 2) + a / (H + alpha * M) ^ 2 * (derH + alpha * der(M)));
    der(Man) = derMan;
    if der(H) > 0 then
        delta = 1;
    else
        delta = -1;
    end if;
    der(M) = ((1 - c) * der(H) * (Man - M)) / ((1 - c) * delta * k - alpha * (Man - M)) + c * der(Man);
    B = 0.000001256* (H + M);
end Hysteretic_curve;

It´s worked now, but results aren´t good. Problem is still with Man. I´ve tried to make it in Simulink and it works fine , but I need it in Modelica ....

There are 0 guests and 0 other users also viewing this topic
You are here: