Signal Transformer Characterization for SPICE
Transformers are darned hard to model for many reasons:
As you have stated, core materials do not have a linear permeability, but vary as they become saturated and also vary with transformer temperature.  If you look at a standard BH curve like the one below you will see the issues.  Each time you magnetize the core in one direction, it retains some magnetism and when you reverse current, energy is required to first demagnetize and then re-magnetize in the opposite direction (called hysteresis).  Since the apparent inductance is related to the slope of the curve, it is dependent on not only on what the instantaneous voltage and di/dt are but what path you took to get there.  Add to this the troubling fact that any DC component of your signal pushes you up or down the curve to a place where the slope of the line is different.  Throw in the resistance of the coils and the capacitance of the wires lying next to each other, and the fact that the hysteresis losses and eddy currents actually get hot and dissipate power (making them non-linear resistive components in your model), and you can understand why attempts at models all have loss factors and other "lumped" type parameters in them.  I guess this is not really an answer, but it wouldn't fit as a comment.