I don't dispute that. When working with L, it's correct to use k which is the geometric mean of k ij and k ji, and when working with circuits one uses L. Flux linkages are certainly not an abstraction to a transformer, though they may well be to a circuit designer. When solving the matrix equations for L's, the K ji's and K ij's get cross-multiplied so once in the realm of L it doesn't matter if they're alike or different since they've disappeared as entities.
I didn't find much help in textbooks or in the literature when I was developing my theory some years ago. Seems like it must be somewhere, but I never found it. I've never published it and don't plan to. I can assure you that it works, y'all can decide whether or not you care to believe that. There was a TDK appnote at the time (1995 or so) that stated something like "rod core transformers cannot be designed analytically, they must be developed empirically." Wrong! I wish I still had a copy of that app note! Maybe using FEMM is empirical, but the results from it wouldn't be of any use without the theory.
Let's try a heuristic argument. By reciprocity, the permance of the flux path linking two windings is the same regardless of which winding is excited. But the permeance of the total flux paths for the windings may be quite different, particularly in loosely coupled situations if the windings are geometrically rather different relative to one another and/or in how they're disposed on a core that doesn't comprise a closed path. Therefore, the ratio of flux linking a pair of windings to total flux for a given winding (k ij) can be quite different. FEMM finds these flux linkages as space integrals so flux partially linking a winding (not the whole winding) gets "partial credit".
Let's look at a simple familiar model. This one is limited to 2 windings while the general case is n windings, but it may help clarify. Consider the classic three-terminal model of a transformer:
------ Z1 - M -----+------ Z2-M------- | M |
--------------------------------------------
L1 is primary inductance, L2 is secondary inductance, M is mutual inductance. In conventional mesh notation, L11 = L1, L22 = L2, L21 = L12 = M. Mr. Guillemin smiles. But the ratios of mutual to self, call them k21 and k12, are M/L1 from the left and M/ L2 from the right. Even though L21 = L12 = M, the ratios are not necessarily the same. From any text, M = k * sqrt (L1 * L2). By simple algebra, k21 = k sqrt (L2/L1) and k12 = k sqrt (L1/L2). The geometric mean, sqrt (K21*K12) = k. QED.
If one is working with an existing transformer, it's definitely easiest to measure and think in terms of inductances. If one is designing loosely-coupled transformers, it's easier to use FEMM to see what flux linkages occur with a given geometry. It's vert easy to get k ji and k ij from a FEMM simulation. One can examine variants that way a lot faster than winding, vacuum-potting and testing actual specimens -- and at 30KV in a small xfmr vacuum-potting is not optional though somewhat messy. Once the flux linkage data is converted to inductances in a math model, Lij = Lji in every case as Guillemin and others assert.
Both approaches are correct, and both work. Transformers I've built accurately to a given geometry usually exhibited circuit parameters within 5% or so of what I expected from FEMM modelling and math modelling. The math models were used to predict the effects of distributed capacitances and hence self-resonances, and they did account for winding resistances. Fringing and all that are non-issues because FEMM analyzes and integrates fluxes as they are determined to actually exist. The resulting values of L ij and k were then used in SPICE to predict circuit behavior with associated electronics. Once actual circuits and transformers were made, they performed at the bench as expected.
The application here was very low cost 30KV 30-watt high-frequency oil-ignition transformers operating at resonance with a capacitor in a self-excited power oscillator. Loose coupling is useful here because it enables the desired high output impedance without requiring additional inductances.
The circuit designer could care less about k ij and k ji, since all he can see and all he need care about is inductances. For closely-coupled designs as used in SMPS, the difference between k ji and k ij must be miniscule in any case since k is "very close to 1". I now know, thanks to the experience and helpful post of Mr Woodgate, that values of .998 are quite achievable so I guess I need not worry about k or k ij or such claptrap for my little SMPS project. I'll just wind the suckuh and snub as necessary. I apologize for the diversion.